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Summary 

This  report  summarizes  the  main  accomplishments  of  a  three-year  research  program  supported 
under  AFOSR  Grant  No.  F49620-93- 1-0177.  The  main  objective  of  this  program  was  to 
investigate  active  suppression  of  detrimental  combustion  instabilities  in  chemical  rockets  by  a 
controlled,  secondary,  combustion  process.  The  program  consisted  of  parallel  theoretical  and 
experimental  efforts;  the  former  developed  the  theoretical  foundation  of  the  investigated  control 
approach  and  the  latter  developed  a  small  scale,  actively  controlled,  gas  rocket  setup  that  was  used 
to  guide  the  development  of  the  investigated  active  control  system  and  demonstrate  its 
effectiveness. 

The  developed  active  control  system  (ACS)  consists  of  a  pressure  transducer  that  continuously 
measures  the  combustor's  pressure,  an  observer  that  analyzes  the  measured  pressure  and 
determines  the  amplitudes,  phases  and  frequencies  of  the  unstable  combustor  modes  in  real  time,  a 
controller  that  provides  each  identified  mode  (by  the  observer)  with  a  phase  shift  and  a  gain  and 
generates  a  control  signal  that  is  sent  to  a  fuel  injector  actuator  that  modulates  the  injection  rate  of  a 
secondary  fuel  stream  into  the  combustor.  This  control  system  is  based  upon  Rayleigh’s  criterion 
and  designed  to  produce  a  secondary,  oscillatory,  combustion  process  within  the  combustor  that  is 
out  of  phase  with  the  combustor  pressure  oscillations,  thus  resulting  in  their  attenuation. 

The  theoretical  efforts  investigated  the  performance  of  an  actively  controlled  rocket  motor  that 
is  prone  to  axial  instabilities.  Two  models  that  investigated  the  control  of  axial  instabilities  were 
developed.  The  first  used  an  approximate,  Galerkin  type,  approach  to  study  the  control  of  linear 
instabilities  and  the  second  used  a  heuristic  model  and  a  numerical  solution  approach  to  investigate 
the  performance  of  an  actively  controlled  rocket.  The  latter  used  a  phenomenological  model  to 
describe  combustor  mixing  processes  and  global,  Arrhenius  type,  global  kinetics  to  describe  the 
combustion  processes.  These  models  were  used  to  predict  the  conditions  within  the  combustor 
with  the  ACS  on  and  off.  The  control  of  both  linear  and  nonlinear  instabilities  was  investigated. 
These  studies  demonstrated  that  the  developed  ACS  can  identify  the  characteristics  of  severe 
instability  in  practically  real  time  and  effectively  damp  each  unstable  mode.  Furthermore, 
numerical  predictions  of  the  combustor  response  to  open  loop  control  were  found  to  be  in  good 
agreement  with  experimental  results.  Finally,  these  studies  developed  improved  means  for  the 
numerical  representation  of  the  injector  face  boundary  condition  and  the  handling  of  various 
sources  (e.g.,  heat  and  mass  addition,  area  change  and  friction)  in  Roe's  Rieman  scheme,  which 
was  used  to  numerically  solve  the  model  equations. 
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The  experimental  efforts  developed  an  actively  controlled  gas  rocket  that  was  subsequently 
used  to  investigate  the  performance  of  the  ACS  in  open  and  closed  loop  control  modes.  The  open 
loop  studies  developed  two  approaches  for  determining  the  frequency  dependence  of  the  fuel 
injector  actuator-combustor  system  response.  The  results  of  these  studies  were  stored  in  the  ACS' 
controller  and  used  to  determine  the  phase  shift  and  gain  that  controller  must  add  to  the  unstable 
combustor  modes  in  closed  loop  control.  Subsequent,  closed  loop,  control  studies  showed  that  the 
investigated  ACS  can  effectively  damp  large  amplitude,  highly  nonlinear,  instabilities  in  periods  of 
the  order  of  40  milliseconds,  which  are  considerably  shorter  than  the  times  reported  by  other 
investigators. 

The  main  contributions  of  this  study  are: 

1 .  the  development  of  an  active  control  approach  based  upon  Rayleigh's  criterion, 

2.  the  development  and  demonstration  of  an  observer  that  determines  the  amplitudes,  phases 
and  frequencies  of  a  prespecified  number  of  unstable,  combustor  modes  in  real  time, 

3.  the  development  of  theoretical  models  for  investigating  active  control  of  linear  and 
nonlinear  instabilities  in  rocket  motors, 

4.  the  development  of  a  fuel  injector  actuator  that  utilizes  a  magnetostrictive  material  to 
modulate  a  secondary  fuel  injection  rate  over  a  0-1,200  Hz.  range,  which  is  wider  than  that  of  any 
known  injector, 

5.  the  demonstration  that  the  developed  fuel  injector  actuator  can  excite  significant  reaction  rate 
heat  release  oscillations  within  the  combustor, 

6.  the  development  of  two  different  approaches  for  determining  the  frequency  dependence  of 
the  response  of  the  secondary  combustion  process  generated  by  the  fuel  injector  actuator,  and 

7.  the  demonstration  that  closed  loop  application  of  the  investigated  ACS  significantly  (e.g., 
by  26  dB.)  and  rapidly  (e.g.,  within  40  milliseconds)  damps  a  rocket  motor  instability  without 
destabilizing  any  stable  modes. 

It  is  believed  that  these  findings  provide  a  foundation  for  guiding  the  development  of  ACS  for 
unstable  rockets  and  other  combustors.  The  results  of  this  program  have  been  demonstrated  to  a 
number  of  companies.  A  major  gas  turbine  manufacturer  is  currently  working  on  adapting  the 
developed  ACS  technology  for  application  in  unstable  gas  turbines,  and  a  fuel  injector  actuator 
based  upon  the  one  developed  under  this  program  is  being  investigated  by  the  Navy  for  application 
in  active  control  of  shipboard  incinerators. 
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Introduction 


This  report  describes  the  results  of  a  three-year  investigation  of  active  control  of  combustion 
instabilities  in  chemical  rockets.  Such  instabilities  are  generally  driven  by  a  feedback-type 
interaction  between  flow  and  combustion  process  oscillations;  energy  supplied  by  an  oscillatory 
combustion  process  produces  oscillatory  heat  release  that  excites  one  or  more  natural  acoustic 
modes  of  the  combustor  whose  oscillations  are  responsible  for  the  periodicity  of  the  combustion 
process.  The  condition  for  driving  an  instability  in  a  combustor  can  be  expressed  by  the  following 
modified  form  of  Rayleigh's  criterion: 

lJrP'(x,t)Q'(x,t)dtdV  ;>{/•  2L,(x,/)<MV 

*  (1) 

where  p',  Q',  Lj,  t,  x,  T  and  V  are  the  combustor  pressure  oscillation,  heat  addition  oscillation,  i-th 
loss  process  (e.g.,  viscous  dissipation,  acoustic  energy  transmission  through  the  nozzle),  time, 
location  within  the  combustor,  period  of  the  oscillation  and  combustor  volume,  respectively.  The 
integral  on  the  left  and  right  hand  side  of  Eq.  1  describes  the  total  driving  and  damping  experienced 
by  the  combustor  oscillations,  respectively.  An  instability  occurs  when  the  inequality  in  Eq.  1  is 
satisfied;  that  is,  the  overall  driving  within  the  combustor  is  larger  than  the  overall  combustor 
damping.  It  is  noteworthy  that  driving  occurs  at  a  given  location  x  when  the  time  integral  on  the 
left  side  of  Eq.  1  is  positive  at  that  location.  This  condition  is  satisfied  when  magnitude  of  the 
phase  difference  <j>  between  the  pressure  and  heat  addition  oscillations  is  less  than  90  degrees.  It 
follows  that  at  locations  where  the  magnitude  of  this  phase  difference  is  larger  than  90  degrees,  the 
oscillatory  heat  addition  process  damps  the  oscillations. 

An  examination  of  Eq.  1  suggests  that  instabilities  could  be  prevented  by  decreasing  the 
magnitude  of  the  integral  on  the  left  and/or  increasing  the  magnitudes  of  various  loss  terms  on  the 
right.  The  former  can  be  attained  by  modifying  the  characteristics  of  the  combustion  process 
and/or  its  interaction  with  the  flow  oscillations.  The  characteristics  of  the  combustion  process 
could  be  changed  by,  for  example,  modifying  the  propellants'  feed  system,  fuel  injectors  and 
propellants'  composition.  On  the  other  hand,  the  loss  terms  on  the  right  hand  side  of  the  Eq.  1 
could  be  increased  by,  for  example,  changing  the  characteristics  of  the  nozzle  and/or  adding 
mechanical  elements  (e.g.,  an  acoustic  liner)  that  dissipate  acoustic  energy  in  the  combustor. 
These  control  approaches  are  generally  referred  to  as  passive,  and  their  effective  implementation 
requires  understanding  of  the  mechanisms  that  drive  the  instability,  the  characteristics  of  the  excited 
oscillations,  and  the  losses  produced  by  various  system  elements. 
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Unfortunately,  passive  control  approaches  have  generally  not  been  satisfactory.  Lack  of 
adequate  understanding  of  the  fundamental  processes  that  control  the  instability  resulted  in 
"solutions"  there  were  only  applicable  to  a  specific  combustor  design  over  a  limited  range  of 
operating  conditions.  Consequently,  different  passive  control  approaches  had  to  be  developed  for 
different  unstable  combustors.  Since  these  passive  control  approaches  were  generally  developed  in 
costly  and  lengthy  trial-and-error  development  programs,  it  became  apparent  that  new  approaches 
for  controlling  combustion  instabilities  are  needed. 

The  considerable  progress  in  the  areas  of  computers,  electronics,  sensors,  actuators  and  control 
theory  in  recent  years  has  renewed  interest  in  the  development  of  active  control  systems  (ACS)  for 
preventing  combustion  instabilities.  An  ACS  can  prevent  the  onset  and/or  damp  an  instability  by 
one  or  more  of  the  following  actions: 

1 .  production  of  heat  addition  oscillations  in  the  combustor  that  are  180  degrees  out  of 
phase  relative  to  the  unstable  pressure  oscillations,  thus  resulting  in  their  attenuation 
(based  upon  Rayleigh's  criterion;  see  Eq.  1  above),  and/or 

2.  interference  with  the  mechanisms  that  drive  the  instability  in  a  way  that  reduces  their 
driving  effectiveness,  and/or 

3.  modification  of  the  system's  boundary  condition(s)  in  a  way  that  produces  one  or  all  the 
following:  (i)  increase  in  the  system’s  damping,  (ii)  modification  of  the  modes  that  can  be 
excited  within  the  system,  and  (iii)  destructive  interference  with  the  mechanisms  that 
drive  the  instability. 

A  typical  ACS  generally  consists  of  a  sensor,  an  observer,  a  controller,  and  an  actuator,  see 
Fig.  1.  The  sensor  (e.g.,  a  pressure  transducer  or  a  photo-multiplier)  continuously  "senses"  the 
conditions  inside  the  unstable  combustor.  The  signal  measured  by  the  sensor  is  transmitted  to  the 
observer  that  determines  the  state  of  the  system.  This  information  is  sent  to  the  controller  where  it 
is  modified,  using  a  specific  control  approach,  and  sent  to  the  actuator  that  "perturbs"  the 
combustor  in  a  controlled  manner.  The  latter  prevents  the  onset  or  damps  the  instability  via  one  or 
more  of  the  actions  described  in  Items  1-3  above. 

The  main  advantages  of  a  given  ACS  are: 

1 .  it  can  prevent  the  onset  and/or  rapidly  attenuate  combustion  instabilities, 

2.  it  can  rapidly  respond  to  changes  in  the  combustor  operating  conditions  and,  thus,  the 
characteristics  of  the  instability,  and 
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3.  it  is  expected  to  be  applicable  (perhaps  with  some  modifications)  to  different  unstable 
combustors,  thus  providing  a  capability  for  effective  control  of  instabilities  in  a  variety  of 
combustion  systems. 

While  prior  investigations  of  active  control  of  combustion  instability  had  demonstrated  the 
considerable  promise  of  this  approach,  a  close  examination  of  the  results  of  these  studies  revealed 
that  much  remained  to  be  done  in  several  areas  before  active  control  can  be  implemented  in  practical 
combustors.  For  example,  past  investigations  generally  filtered  the  measured  signal  before 
sending  it  to  the  controller  where  it  was  amplified  and  phase  shifted  before  being  sent  to  the 
actuator.  The  use  of  a  filter  requires,  however,  apriori  knowledge  of  the  frequency  of  the 
instability.  While  this  frequency  can  be  readily  determined  in  advance  in  experimental  combustors, 
it  is  generally  not  known  apriori  in  practical  combustors,  and  may  vary  in  response  to  changes  in 
combustor  operating  conditions.  Consequently,  a  practical  ACS  should  possess  capabilities  for 
determining  in  real  time  the  frequencies  of  the  unstable  combustor  modes.  Another  problem  is  the 
shortcomings  of  actuators  such  as  loudspeakers,  mechanical  valves  and  fuel  injectors  that  had  been 
utilized  by  other  investigators.  For  example,  loudspeakers  generally  don't  have  the  power  required 
for  stabilizing  practical  combustors,  capabilities  for  continuous  operation  without  failure,  and 
"hardware"  that  can  survive  in  hostile  combustor  environments.  Similarly,  mechanical  valves  are 
complex,  heavy  and  limited  to  low  frequency  applications.  It,  thus,  became  apparent  at  the  outset 
of  this  program  that  a  fuel  injector  that  can  supply  an  oscillatory  fuel  flow  rate  into  the  combustor 
to  drive  pressure  oscillations  at  a  desired  frequency  and  phase  and/or  modify  the  primary 
combustion  process  in  a  manner  that  reduces  its  driving  effectiveness  should  be  used  as  the 
actuator  in  an  ACS  for  damping  combustion  instabilities.  The  automotive  fuel  injectors  that  were 
utilized  in  past  studies  are  limited,  however,  to  low  frequencies  and  low  fuel  flow  rates  and, 
consequently,  are  not  suitable  for  use  in  practical  ACS. 

In  an  effort  to  advance  the  state  of  the  art  and  develop  capabilities  that  would  lead  to  the 
development  of  practical  active  control  systems  for  unstable  rocket  (and  other)  combustors,  efforts 
under  this  program  focused  on  developing  the  following  capabilities: 

1 .  an  observer  that  can  analyze  the  sensor's  signal  and  determine  in  the  amplitudes, 
frequencies  and  phases  of  up  to  five  unstable  combustor  modes  in  virtually  real  time, 

2.  a  fuel  injector  actuator  that  can  modulate  the  flow  rate  of  a  gaseous  fuel  with  significant 
amplitudes  over  a  0-1,200  Hz.  frequency  range,  whose  upper  limit  is  considerably 
higher  than  the  maximum  frequency  of  any  known  (e.g.,  automotive)  fuel  injector. 
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3 .  a  small  scale,  actively  controlled,  variable  length,  gas  rocket  that  exhibits  axial 
instabilities,  similar  to  those  observed  in  full  scale  rocket  motors,  over  a  wide  range  of 
frequencies  (e.g.,  100-1800  Hz.),  and 

4.  a  model  of  an  actively  controlled  gas  rocket  that  can  be  used  to  theoretically  investigate 
the  performance  of  various  ACS  and  control  strategies. 

The  research  activities  and  accomplishments  of  this  program  are  described  in  the  remainder  of 
this  report.  These  are  described  in  the  following  order:  1.  the  principles  of  operation  of  the 
investigated  ACS,  2.  theoretical  investigations  of  the  performance  of  unstable  rocket  motors 
actively  controlled  with  the  investigated  ACS,  3.  the  developed  experimental  setup,  4.  open  loop 
active  control  performance  of  the  investigated  ACS,  5.  closed  loop  active  control  of  combustion 
instabilities  with  the  investigated  ACS,  and  6.  summary  and  conclusions  of  the  results  of  this 
program. 
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The  Investigated  Active  Control  System 

At  the  onset  of  this  program  it  was  recognized  that  an  ACS  capable  of  preventing  detrimental 
rocket  motor  instabilities  will  require  capabilities  for: 

1.  identification  of  the  characteristics  of  all  the  unstable  combustor  modes  in  virtually  real  time 
(i.e.,  before  the  instability  results  in  rocket  malfunction)  using  limited  sensor  data  about  the 
combustor  performance, 

2.  determination  of  the  appropriate  phase  and  gain  that  must  be  "added"  by  the  controller  to 
each  unstable  mode  and  generation  of  an  appropriate  control  signal  for  the  actuator,  and 

3.  means  for  introducing  controlled  "perturbations"  into  the  combustor  that  will  attenuate  the 
unstable  combustor  modes. 

To  provide  these  capabilities,  the  developed  ACS  uses  a  pressure  transducer  installed  near  the 
combustor's  injector  as  its  sensor,  because  the  anti-nodes  of  most  potentially  unstable  axial 
pressure  oscillations  are  likely  to  occur  near  this  location.  This  maximizes  the  likelihood  that  the 
measured  pressure  will  include  contributions  from  all  unstable  combustor  modes.  The  measured 
pressure  is  then  transmitted  to  an  observer  whose  task  is  to  determine,  in  virtually  real  time,  the 
amplitudes,  frequencies  and  phases  of  all  combustor  modes  that  significantly  contribute  to  the 
instability.  The  identified  characteristics  of  the  unstable  modes  are  transmitted  to  a  controller  where 
each  mode  is  provided  with  a  gain  and  a  phase  shift  that  will  result  in  the  attenuation  of  the 
corresponding  mode  within  the  combustor.  The  controller  combines  the  modified  modes  into  a 
single  signal  that  is  sent  to  a  fuel  injector  actuator  that  can  modulate  the  injection  rate  of  a  gaseous 
fuel  into  the  combustor  over  a  0-1,200  Hz.  frequency  range.  The  fuel  injection  modulations 
produce  reaction  rate  and  heat  addition  oscillations  within  the  combustor  at  the  frequency  of  the 
fuel  flow  rate  modulation.  These  heat  addition  oscillations  should,  according  to  Rayleigh’s 
criterion,  damp  the  instability  if  the  magnitude  of  the  phase  difference  between  these  heat  addition 
and  combustor  pressure  oscillations  is  smaller  than  ninety  degrees. 

The  Observer 

Since  the  developed  observer  plays  a  key  role  in  the  investigated  ACS  and  it  is  repeatedly 
referred  to  in  the  remainder  of  this  report,  its  principles  of  operation  are  presented  in  this  section. 

The  observer  analyzes  the  measured  pressure  to  determine  in  real  time  the  amplitudes,  phases, 
and  frequencies  of  the  unstable  combustor  modes.  It  assumes  that  the  measured  combustor 
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pressure  p(t)  can  be  expressed  as  the  following  sum  of  combustor  modes  that  may  not  be 
harmonics  of  one  another  and  may  have  frequencies,  amplitudes  and  phases  that  slowly  vary  with 
time. 


n=k 


p(t)=  1LSnsm(cont)  +  Cncos(o)nt ) 

n= 1 


(2) 


The  observer  determines  the  characteristics  of  the  dominant  (i.e.,  largest  amplitude)  mode  by 
solving  the  integrals 

2  1  2  1 

^n(0  =  —  Jsin {cont)p{t)dt\  Cn  (t)  =  — —  jcos(cont)p(t)dt  (3) 

1  n  t~Tn  1  n  t-Tn 

where  the  oscillation  period  Tn  and  frequency  COn  are  related  by  Tn-2nl  COn  and  are  not 

known  apriori.  It  can  be  shown  that  the  above  integrals  can  replaced  by  following  recursive 
formulae: 

Sn(t  +  At)  =  Sn(t)  +  -^-(p(t  +  At)~  p(t-Tn  +  At))  •  sin(cont)  ■  At  (4-a) 

2 

Cn(t  +  At)=  Cn(t)  +  — ~(p(t  +  At)-  p(t-Tn  +  At))  ■  cos(cont)  ■  At  (4-b) 

n 


whose  solution  requires  little  computational  effort1 

The  unknowns  Sn(t),  Cn(t)  and  COn  are  determined  in  a  rapidly  converging  iterative  solution 
procedure.  Initially,  a  value  for  COn  (and,  thus,  Tn)  is  assumed  and  substituted  into  Eqs.  4,  which 
are  then  solved  for  the  "corrected"  coefficients  of  Sn(t+Dt)  and  Cn(t+Dt).  The  calculated 

coefficients  are  then  substituted  into  the  following  relationship1  that  determines  a  corrected  value  of 
the  frequency: 


[M2] 


t+At 


CO 


f  (Js  dc  *\ 

co  —  sin(coO  +— cos(cot) 
\dt _ at _ 

S  •  cos(cot) -  C •  sin(atf) 


it 


(5) 


where  the  subscript  'n'  has  been  omitted  for  simplicity.  The  values  of  CO  on  the  left  and  right  sides 
of  Eq.  5  are  the  "corrected"  and  "previous  calculated"  values  of  the  frequency,  respectively.  Once 
a  corrected  value  of  CO  is  obtained  from  Eq.  5,  the  corresponding  period  T  is  substituted  into  Eqs. 


9 


4  to  obtain  improved  values  of  S(t)  and  C(t).  This  procedure  rapidly  converges  into  the  "final" 
values  of  the  unknowns  CO,  S  and  C,  which  are  then  used  to  determine  he  amplitude,  phase  and 
frequency  of  the  dominant  mode.  This  procedure  can  also  track  "slow"  variations  of  these 
quantities.  Once  the  characteristics  of  the  dominant  mode  are  known,  the  expression  describing  its 
time  dependence  is  subtracted  from  the  measured  pressure  p(t)  and  the  above  procedure  is  repeated 
to  determine  the  characteristics  of  the  "next"  dominant  mode  within  the  remaining  signal. 

The  above  described  procedure  can  be  repeated  to  determine  the  characteristics  of  as  many 
modes  as  desired.  Clearly,  the  need  to  attain  real  time  identification  of  the  characteristics  of  an 
instability  would  limit  the  number  (e.g.,  two  or  three)  of  modes  that  can  be  identified.  This  is  not, 
however,  a  problem  as  most  instabilities  are  generally  dominated  by  one  or  two  modes. 

Figures  2-a  through  2-e  demonstrate  the  ability  of  the  observer  to  identify  in  real  time  the 
hierarchy  of  modes  in  an  unstable  combustor.  The  observer  was  provided  with  the  pressure  signal 
shown  in  Fig.  2-a,  which  was  measured  in  the  gas  rocket  motor  that  was  developed  under  this 
program  when  it  experienced  combustion  instability.  Figure  2-a  shows  that  the  instability 
"switches"  from  low  to  high  frequency  oscillations  between  .05  and  .064  seconds.  This  signal 
was  analyzed  by  the  observer  and  the  computed  frequencies  and  amplitudes  were  used  to  determine 
the  time  dependence  of  the  two  most  dominant  combustor  modes,  see  Figs.  2-b  through  2-e. 
Figures  2-b  shows  that  the  frequency  of  the  nearly  sinusoidal,  dominant,  mode  abruptly  increases 
around  .065  seconds  to  that  of  the  second  mode,  while  Fig.  2-c  shows  that  the  second  mode 
oscillates  with  a  higher  frequency  and  its  amplitudes  decreases  to  practically  zero  around  .065 
seconds.  Figures  2-b,c  also  show  that  the  observer  can  simultaneously  track  the  behavior  of  both 
modes  in  real  time.  Figure  2-d  presents  the  calculated  time  dependence  of  the  frequencies  of  the 
two  observed  modes.  It  shows  that  the  observed  frequency  of  the  dominant  mode  changed  from 
650  to  1250  Hz.  within  only  three  milliseconds.  Finally,  Fig.  2-e  compares  the  measured  pressure 
p(t)  with  that  obtained  by  synthesis  of  the  two  observed  modes.  It  clearly  shows  that  the  observer 
can  "reproduce"  the  input  pressure  in  real  time  and  with  high  fidelity. 

During  the  last  year  of  this  program  a  theoretical  study  aimed  to  extend  the  initially  developed 
observer  approach  to  provide  a  capability  for  simultaneous  identification  of  several  unstable  modes 
was  undertaken.  This  study  is  presented  in  Appendix  A.  While  the  current  observer  determines 
the  characteristics  of  the  unstable  modes  in  a  hierarchical  manner  (i.e.,  the  most  dominant  mode  is 
identified  first,  then  the  next  mode  dominant  mode  is  identified  and  so  on),  the  new  approach 
simultaneously  identifies  the  characteristics  of  several  unstable  combustor  modes.  It  is  expected 
that  the  results  of  this  study  will  increase  the  accuracy  and  robustness  of  the  observer  and  increase 
the  likelihood  that  its  output  will  converge  to  the  correct  results. 
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Theoretical  Studies 


The  performance  of  the  investigated  ACS,  see  Fig.  1,  was  initially  theoretically  studied  by 
solving  the  system  of  one  dimensional,  unsteady,  conservation  equations  that  modeled  the  flow  in 
an  unstable  combustor^.  The  model  assumed  that  the  instability  is  driven  by  a  simple  linear 
feedback  between  the  pressure  and  combustion  process  heat  addition  oscillations.  The  unstable 
combustor  model  was  then  "equipped"  with  the  developed  ACS  and  used  to  theoretically 
investigate  the  effectiveness  of  the  developed  ACS.  Numerical  simulations  showed  that  the 
developed  ACS  can  rapidly  damp  a  variety  of  instabilities  including,  for  example,  combustor 
acoustic  mode  oscillations  with  amplitudes  equaling  twenty-five  percent  of  the  mean  combustor 
pressure,  which  were  damped  within  a  few  milliseconds  after  activating  the  ACS. 

Subsequent  studies  developed  a  heuristic  model  of  an  unstable  gaseous  rocket  combustor  and 
an  improved  numerical  approach  for  solving  the  model  equations.  The  conservation  equations 
solved  by  this  model  are  summarized  in  Fig.  3.  Since  the  model  assumes  that  the  combustor  flow 
is  one  dimensional  and  inviscid,  and,  thus,  cannot  account  for  turbulent  mixing,  a 
phenomenological  mixing  model  that  describes  the  mixing  of  the  premixed  reactants  with  the  hot 
combustion  products  has  been  incorporated  into  the  conservation  equations  (i.e.,  see  terms 
proportional  to  l/Xmi*).  This  model  also  accounts  for  the  mixing  of  the  secondary  fuel,  injected  by 
the  fuel  injector  actuator,  with  the  combustor  flow  before  its  reaction.  The  magnitude  of  the  spatial 
dependence  of  the  mixing  process  is  controlled  by  the  function  W(x),  see  Fig.  3,  which  depends 
upon  a  prespecified  mixing  length  lmix-  Once  a  combustible  mixture  is  formed,  the  combustion 
process  heat  release  rate  is  controlled  by  global,  Arrhenius  type,  methane  kinetics,  which  is  the  fuel 
used  in  the  experimental  phase  of  the  program. 

Once  it  had  been  shown  that  the  heuristic  model  predicts  unstable  rocket  operation  of  ranges  of 
design  and  operating  parameters,  the  model  was  "equipped"  with  the  investigated  ACS  and  the 
resulting  model  was  used  to  investigate  the  behavior  of  an  actively  controlled  rocket  motor. 
Typical  results  obtained  in  this  study  are  presented  in  Figs.  4  and  5.  Figure  4  shows  predictions  of 
the  motor's  performance  under  open  loop  control  excitation.  It  shows  predicted  time  dependence 
of  pressure  and  heat  release  oscillations  within  the  rocket  combustor  when  the  fuel  injector  actuator 
modulated  the  secondary  fuel  injection  rate  at  specific  frequency.  The  results  show  that  the 
actuator  can  excite  significant  pressure  and  heat  release  oscillations  within  the  combustor, 
indicating  that  it  could  serve  as  an  actuator  in  closed  loop  control  of  rocket  motor  combustion 
instabilities.  An  analysis  of  the  time  dependence  of  the  results  shows  that  the  heat  release 
oscillations  lead  the  pressure  oscillations  and  that  the  magnitude  of  the  phase  difference  is  always 
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smaller  than  ninety  degrees,  in  agreement  with  Rayleigh's  criterion.  This  phase  difference 
decreases  as  the  amplitude  of  the  excited  oscillations  increases,  and  it  approaches  ninety  degrees  as 
the  growth  grate  of  the  oscillations  decreases,  indicating  that  less  driving  is  needed  when  the 
amplitude  growth  rate  goes  to  zero. 

Figure  5  presents  predictions  of  the  effect  of  closed  loop  control  upon  the  oscillations  in  an 
unstable  rocket  combustor.  The  two  top  plots  show  the  onset,  growth  and  decay  of  the  instability 
before  and  after  the  activation  of  the  ACS  at  t~  03  seconds.  These  plots  show  that  the  ACS 
significantly  attenuated  the  instability.  The  middle  plot  compares  the  predicted  and  "observed"  (by 
the  ACS'  observer)  pressure  oscillations.  It  shows  that  the  developed  observer  can  indeed 
"identify"  the  characteristics  of  the  combustor  pressure  oscillations  in  practically  real  time.  The 
second  plot  in  the  middle  demonstrates  the  observer's  ability  to  "follow"  in  real  time  the  variations 
in  the  frequency  of  the  instability.  The  plot  on  the  bottom  left  presents  an  "amplified"  view  of  the 
oscillations  shown  in  the  top  left  plot  after  activating  the  ACS.  It  shows  that  the  amplitude  of  the 
damped  combustor  oscillations  is  negligible.  Finally,  the  plot  on  the  bottom  right  shows  the  time 
dependence  of  the  control  signal  o  the  actuator.  In  summary,  the  results  presented  in  Fig.  5  predict 
that  the  investigated  ACS  will  effectively  attenuate  rocket  motor  combustion  instabilities. 

In  a  parallel  study,  Roes  Rieman  Solver^  was  modified  for  applications  in  numerical 
simulations  of  unstable  combustors.  Specifically,  Roe’s  Rieman  scheme  was  modified  to  properly 
account  for  the  presence  of  source  terms  (e.g.,  the  combustion  process  heat  addition,  friction)  in 
the  numerical  solution  of  the  conservation  equations.  Furthermore,  numerical  schemes  that 
improve  the  representation  of  the  injector  face  and  nozzle  boundary  conditions  were  developed  for 
incorporation  into  the  numerical  solution  scheme. 

Finally,  the  feasibility  of  applying  open  loop,  nonlinear,  high  frequency,  vibrational  control  to 
damp  instabilities  in  mechanical  systems  was  theoretically  investigated4.  Such  control  has  the 
advantage  that  it  does  not  require  accurate  determination  of  the  state  of  the  controlled  system  and 
the  data  processing  capabilities  required  in  conventional  feedback  control.  This  paper  extends 
previous  linear  investigations  of  this  type  of  control  to  include  nonlinear  vibrational  controllers.  It 
shows  that  a  nonlinear  vibrational  controller  can  stabilize  a  system  even  if  the  Jacobian  matrix  has  a 
positive  trace.  A  copy  of  this  paper,  which  was  submitted  for  publication  in  the  IEEE  Journal,  is 
provided  in  Appendix.B. 

Experimental  Setup 

The  experimental  efforts  led  to  the  development  of  an  actively  controlled,  small  scale,  gas 
rocket  motor  setup  shown  in  Fig.  6.  It  consists  of  a  reactants  feed  system,  a  combustor  section,  a 
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nozzle,  a  secondary  fuel  injection  system  (i.e.,  the  fuel  injector  actuator),  a  main  control  panel,  and 
a  computer-based  controller  of  the  secondary  fuel  actuator.  The  setup  also  includes  transducers 
and  photomultipliers  that  measure  the  combustor  pressure  and  reaction  rate,  respectively,  and  a 
computer-based  data  acquisition  system. 


A  primary  reactants  stream,  consisting  of  premixed  air  and  methane,  enters  the  combustor 
through  the  injector  orifices.  The  air  and  methane  are  separately  introduced  into  a  mixing  port 
through  calibrated,  choked,  orifices  that  provide  means  for  measuring  the  flow  rate  of  each 
reactant.  These  flow  rates  are  manually  controlled  by  setting  the  pressure  upstream  of  each  orifice 
to  a  desired  value.  The  flows  through  the  injector  orifices  are  choked  to  prevent  feedback  between 
the  combustor  and  reactants  supply  lines.  The  injected  reactants  jets  are  oriented  twenty  degrees 

relative  to  the  combustor  axis,  which  forces  the  primary  combustion  process  to  stabilize  near  the 
combustor  walls. 

The  combustor's  exhaust  nozzle  is  choked  at  all  test  conditions.  Under  typical  operating 
conditions,  the  combustor  pressure,  premixed  reactants  line  pressure  and  fuel  and  air  supply  lines 
pressures  are  45, 90  and  200  psi,  respectively.  The  combustor  pressure  was  always  sufficient  for 
choking  the  exhaust  nozzle.  The  supply  pressure  to  the  secondaiy  fuel  injector  is  450  psi.  Under 
this  operating  condition,  the  flow  rates  of  the  primary  stream  of  premixed  fuel  and  oxidizer  and 

secondary  fuel  stream  are  10  and  0.1  gram/sec.,  respectively,  resulting  in  a  total  combustor  power 
output  of  55  kW. 

The  combustor  has  a  large  length  to  diameter  ratio  to  prevent  excitation  of  transverse  modes. 
The  combustor  length  can  be  changed  by  adding  or  removing  pipe  sections.  This  capability 
permits  the  investigation  of  axial  mode  instabilities  having  different  frequencies;  high  frequency 
instabilities  are  excited  in  short  combustors  and  vice  versa.  At  its  shortest  length,  the  combustor's 
fundamental  mode  frequency  is  around  1800  Hz.  while  at  its  longest  configuration  the  fundamental 
mode  frequency  is  200  Hz.  During  operation,  the  combustor  is  immersed  in  a  bath  with  running 
water  to  cool  the  combustor. 

As  shown  in  Fig.  6,  the  axis  of  the  investigated  combustor  is  "broken"  to  provide  a  location  for 
installing  the  window,  optics  and  photomultiplier  required  to  measure  the  total  (integrated) 
combustion  zone  radiation.  This  system  measures  the  total  CC  or  CH  radicals  radiation  from  the 
combustion  region,  which  is  proportional  to  the  combustion  process  heat  release  rate.  A  second, 
side  view,  window,  which  is  installed  near  the  injector,  is  used  for  optical  measurements  and 
visualization  of  the  combustion  zone. 
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The  actively  controlled  secondary  fuel  injector  actuator  introduces  an  oscillatory  fuel  flow  into 
the  combustion  zone.  High  pressure  fuel  is  supplied  to  the  actuator  and  forced  through  an  annular 
orifice  between  the  outer  wall  of  a  needle's  wider  diameter  base  and  its  seat.  A  magnetostrictive 
actuator  is  attached  to  the  needle  and  used  to  change  the  needle's  base  position  relative  to  its  seat  in 
response  to  changes  in  an  electric  control  signal.  The  acoustic  impedance  and  pneumatic  resistance 
of  the  actuator  s  cavities  that  carry  the  fuel  from  the  supply  line  to  the  combustor  were  sized  to 
maximize  the  fuel  flow  rate  oscillations  at  the  injector’s  exit  and  minimize  the  effect  of  the 
combustor  pressure  oscillations  upon  the  actuator's  performance  over  a  wide  frequency  range.  As 
a  result,  this  actuator  can  modulate  the  fuel  flow  rate  over  a  0-1, 200  Hz.  frequency  range  when 
pressure  oscillations  are  present  in  the  combustor  without  a  significant  attenuation  or  a  phase  delay. 
Furthermore,  this  actuator  can  produce  fuel  flow  rate  amplitudes  of  the  order  of  0.2  grams/sec., 
which,  in  principle,  can  produce  20  kW  peak  to  peak  heat  release  rate  oscillations. 

The  electric  signal  to  the  actuator  consists  of  a  steady  and  an  oscillating  component,  which 
control  the  magnitudes  of  the  steady  and  oscillating  flow  rates  through  the  fuel  injector  actuator, 
respectively.  Both  signals  are  generated  in  the  control  computer,  separately  amplified,  and  then 
combined  into  a  single  control  signal  that  is  fed  to  the  actuator. 

Two  different  actuators  and  injection  systems  were  developed  and  tested  during  this  study. 
Each  injection  system,  see  Figs.  7-a,b,  supplied  the  combustor  with  a  primary  stream  of  premixed 
reactants  and  a  secondary  fuel  stream.  The  latter  was  used  to  control  the  instability  by  generating  a 
secondary,  oscillatory,  combustion  process  within  the  combustor.  The  actuators  were  similar  in 
design,  but  the  annular  cross  sectional  area  of  the  second  actuator  was  between  two  to  three  times 
larger  than  that  in  the  first  actuator.  For  further  reference,  the  small  and  large  actuators  will  be 
denoted  as  ACT1  and  ACT2,  respectively,  and  the  injection  systems  shown  in  Figs.  7-a,b  will  be 
denoted  by  INJ1  and  INJ2,  respectively.  The  investigated  injector/actuator  configurations,  whose 
performance  is  described  later  in  this  report,  will  be  described  by  the  notation  INJi/ACTi  where  the 
index  i=l,2  describes  the  injector  and  actuator  used  in  the  study.  It  should  be  pointed  out  that  the 
ACT2  actuator,  which  was  developed  later  in  the  program,  could  be  retrofitted  to  both  injection 
systems  while  the  ACT1  actuator  could  be  only  installed  on  the  INJ1  injector. 

The  second  injector,  INJ2,  and  actuator,  ACT2,  were  developed  and  investigated  in  a 
cooperative  effort  with  a  major  gas  turbine  manufacturer  that  is  interested  in  developing  capabilities 
for  active  control  of  instabilities  in  large  scale  gas  turbines. 
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Open  Loop  Studies 

As  stated  above,  effective  application  of  the  investigated  ACS  in  closed  loop  control  of 
combustion  instabilities  requires  knowledge  of  the  gain  and  phase  shift  that  the  controller  must  add 
to  each  combustor  unstable  mode.  Since  each  unstable  mode  oscillates  at  a  different  frequency  and 
these  frequencies  may  span  a  wide  range,  the  frequency  dependence  of  the  phase  shift  and  gain  that 
the  controller  must  add  to  the  unstable  modes  must  be  known  over  a  wide  frequency  range.  Since, 
unfortunately,  the  characteristics  of  the  secondary  combustion  processes  cannot  be  accurately 
modeled  to  predict  the  frequency  dependence  of  its  phase  shift  and  gain,  these  quantities  must  be 
experimentally  determined  in  open  loop  tests.  The  objective  of  the  open  loop  tests  is  to  measure 
the  phase  delay  between  the  combustion  process  heat  release  and  actuator  control  signal  oscillations 
and  the  ratio  of  the  amplitudes  of  the  combustion  process  heat  release  and  control  signal 
oscillations,  which  is  referred  to  as  the  gain.  Together,  the  measured  phase  shift  and  gain  describe 
the  secondary  combustion  process  response  (SCPR). 

Another  critically  important  reason  for  investigating  the  SCPR  is  to  determine  whether  the 
developed  fuel  injector  actuator  can  excite  combustion  process  heat  release  and  pressure 
oscillations  within  the  combustor  of  sufficient  magnitude  (i.e.,  capable  of  affecting  the  unstable 
combustor  oscillations) .  If  the  developed  fuel  injector  actuator  could  not  excite  heat  addition 
oscillations  of  sufficient  magnitude  within  the  combustor,  it  would  not  be  suitable  for  application  in 
closed  loop  control  of  combustion  instabilities.  In  this  case,  a  different  fuel  injector  (or  another 
type  of  actuator)  will  have  to  be  developed  and/or  investigated. 

Two  approaches  for  determining  the  SCPR  frequency  dependence  were  developed  and 
investigated  under  this  program.  The  first  used  a  very  short  combustor  whose  fundamental, 
longitudinal,  acoustic  mode  frequency  was  around  1800  Hz.,  which  was  significantly  above  the 
upper  limit  of  the  frequency  range  of  the  open  loop  studies.  Consequently,  the  operation  of  the 
combustor  was  "quiet"  over  the  investigated  frequency  range  and  the  magnitudes  of  driven 
disturbances  (by  the  fuel  injector  actuator)  were  significantly  above  the  combustor  noise  level. 
When  tests  were  conducted  in  this  short  combustor,  high  frequency  reaction  rate  and  pressure 
oscillations  driven  by  a  combustor  instability  could  be  readily  observed  in  the  measured  pressure 
and  radiation  data  and  clearly  distinguished  from  the  "single"  frequency  oscillations  driven  by  the 
fuel  injector  actuator.  Using  this  combustor,  the  characteristics  of  the  SCPR  were  determined  from 
the  measured  pressure  and  a  relationship  between  these  pressure  oscillations  and  reaction  rate 

Global  radiation  measurements  could  not  be  used  with  this  combustor  because  the  available  sideview 
window,  see  Fig.  6,  could  not  "view"  the  whole  combustion  zone. 
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oscillations,  which  was  provided  by  a  model  of  the  combustor  oscillations  (which  is  derived 
below).  The  "short"  combustor  also  included  a  side  view  window,  similar  to  the  one  shown  in 
Fiji-  6,  which  was  used  to  measure  CH  radicals  radiation  from  the  combustion  region. 
Unfortunately,  this  window  could  not  "view"  the  whole  combustion  region,  and  the  measured 
radiation  could  only  provide  qualitative  description  of  the  characteristics  of  the  secondary 
combustion  process  oscillations. 

The  following  analysis  describes  the  derivation  of  a  model  that  predicts  the  behavior  of  small 
amplitude  (i.e.,  linear)  oscillations  in  the  short  combustor.  This  model  was  used  to  determine  the 
SCPR  from  measured  pressure  data. 

Neglecting  kinetic  energy  terms,  the  energy  equation  for  a  combustor  with  uniform  properties 
can  be  expressed  in  the  following  form: 

—  J^pedV  =  qcomb  +  maha  +mfhf-  mehe  (6) 


where  the  terms  from  left  to  right  are  the  rate  of  change  of  the  combustor  internal  energy, 
oscillatory  heat  release  supplied  by  the  secondary  combustion  process,  the  air  and  fuel  enthalpy 
fluxes  into  the  combustor  and  the  enthalpy  flux  exiting  through  the  choked  nozzle,  respectively. 
Assuming  a  perfect  gas  behavior  and  that  the  flow  through  the  choked  nozzle  is  quasi  steady,  the 
following  relationships  can  be  derived: 

|WV=(VCv/i?)^;  mehe  =  (KNp / -jT)CpT  (7) 


Substituting  Eqs.  (7)  into  Eq.  (6),  assuming  that  each  dependent  variable  consists  of  a  steady 
component  and  a  small  amplitude  perturbation  (e.g.,  p  =  p  +  p' ),  and  linearizing  the  resulting 
energy  equation  yields  the  following  linear  form  of  the  energy  equation 


=J_  =  l  ,  rdp  Tj_'  VTaP 

mCpT  p  p  dt  2 T'  TClm 


(8) 


where  T  is  a  characteristic  time,  Ca  is  the  speed  of  sound  and  the  subscript  'a'  denotes  the  state  of 
the  incoming  air.  Since  the  temperature  perturbations  generated  within  the  combustion  process  are 
expected  to  decay  before  they  reach  the  nozzle  entrance  at  high  frequencies,  the  term  involving  T'e 
is  neglected  in  Eq.  (8),  yielding  the  following  relationship 
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for  determining  the  heat  release  q  oscillations  from  measured  pressure  oscillations  at  high 
frequencies. 

In  the  second  approach  for  determining  the  SCPR,  the  combustor  shown  in  Fig.  6  was 
employed.  Since  this  combustor  is  considerably  longer  than  the  "short"  combustor  that  was  used 
in  the  first  study,  several  natural  longitudinal  acoustic  modes  of  the  combustor  could  be  excited  in 
the  range  of  investigated  SCPR  frequencies.  For  example,  the  fundamental  mode  frequency  of 
this  combustor  was  370  Hz.  Consequently,  small  heat  addition  oscillations  at  one  of  the  natural 
acoustic  mode  frequencies  of  the  combustor  could  produce  significant  pressure  oscillations  at  this 
frequency.  Consequently,  when  the  fuel  injector  actuator  is  operated  at  a  non-resonant  frequency 
of  this  combustor,  it  excites  significant  secondary  combustion  process  oscillations  within  the 
combustor  but  only  small  amplitude  pressure  oscillations.  Thus,  the  phase  and  amplitude  of  the 
secondary  combustion  process  heat  release  oscillations  can  be  determined  at  non-resonant 
frequencies  of  the  combustor  by  direct  measurements  of  the  characteristics  of  the  total  radiation 
from  the  combustion  region.  These  were  directly  measured  with  a  photomultiplier  through  the 
window  on  the  slanted  combustor  section  in  Fig.  6. 

Figure  8  shows  an  example  of  the  time  dependence  and  spectra  of  pressure  and  CH  radiation 
signals  measured  in  the  combustor  shown  in  Fig.  8  when  the  secondary  fuel  injection  rate 
oscillated  at  610  Hz.  Figure  8-a  shows  that  the  pressure  spectrum  is  dominated  by  spikes 
representing  three  unstable  combustor  modes  at  370,  740  and  1 1 10  Hz.  In  general,  if  a  single 
mode  is  driven  by  an  unstable  combustion  process,  the  amplitudes  of  its  harmonics,  which  are 
excited  by  nonlinear  processes,  are  smaller  than  the  amplitude  of  the  most  unstable  lowest 
frequency  mode.  The  pressure  spectrum  in  Fig.  8-a  indicates  that  this  was  not  the  case  in  this 
experiment  as  the  amplitudes  of  the  higher  frequency  spikes  are  larger  than  the  amplitude  of  the 
fundamental  mode.  This  suggests  that  each  of  the  observed  pressure  oscillations  was 
independently  driven  by  the  combustion  process. 

An  examination  of  the  pressure  spectrum  in  Fig.  8  also  shows  a  small,  but  clearly  visible, 
spike  at  610  Hz.,  the  frequency  of  the  actuator.  In  contrast,  only  one,  large  amplitude,  spike  is 
present  in  the  spectrum  of  the  CH  radiation  signal  at  610  Hz.  The  dominance  of  the  spike  of  the 
driven  heat  release  oscillations  in  the  CH  radiation  spectrum  indicates  that  the  fuel  injector  actuator 
could  readily  excite  heat  release  oscillations  in  this  combustor  that  are  significantly  larger  than  any 
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other  heat  release  oscillation  when  the  frequency  of  the  fuel  flow  rate  modulation  was  not  close  to 
any  natural  acoustic  mode  frequency  of  the  combustor. 

Figure  8  shows  that  FFT  can  be  used  to  identify  a  small  amplitude,  periodic,  pressure  signal  in 
the  neighborhood  of  large  amplitude  spikes.  It  was  found,  however,  that  an  FFT  analysis  cannot 
accurately  determine  the  phase  of  such  a  small  signal.  Unlike  the  FFT  determination  of  the 
amplitude,  the  FFT  determination  of  the  phase  is  highly  sensitive  to  the  frequency  at  which  the 
FFT  integral  is  evaluated.  Since  the  FFT  evaluation  is  performed  at  a  discrete  number  of 
frequencies,  it  may  not  be  performed  at  the  exact  frequency  of  the  secondary  combustion 
oscillations,  thus  increasing  the  likelihood  of  inaccurate  phase  determination. 

To  overcome  the  shortcomings  of  the  FFT,  a  MATLAB-based  software  that  employs  an 
ensemble  averaging  technique  was  developed  to  accurately  determine  the  amplitude  and  phase  of 
the  secondary  heat  release  oscillations.  The  basic  steps  of  this  approach  are: 

1 .  The  input  to  the  actuator,  (usually  a  current)  is  chosen  as  a  reference  signal. 

2 .  Measured  data  points  are  correlated  with  the  reference  signal  by  referencing  the  time  of  each 
data  point  to  a  specific  phase  in  the  period  of  the  reference  signal,  thus  "collapsing"  all  the 
measured  data  points  into  one  period. 

3.  A  "moving"  narrow-width  window  is  used  to  obtain  the  ensemble  average  of  the  time 
dependence  of  the  collected  data  over  the  period  of  the  reference  signal.  It  should  be  noted  that  the 
time  dependence  of  the  calculated  average  does  not  necessarily  assume  a  sinusoidal  shape. 

4)  The  ensemble-averaged  line  is  curve-fitted  with  a  sinusoidal  signal  to  obtain  the 
corresponding  amplitude  and  phase  of  the  measured  data. 

Figure  9  presents  typical  results  obtained  with  this  data  reduction  approach,  using  the  data 
presented  in  Fig.  8.  The  top  to  bottom  plots  on  the  left  side  of  Fig.  9  show  measured  data  points 
describing,  respectively,  the  actuator  control  signal  (current),  actuator  needle  displacement, 
combustor  pressure  and  combustion  process  CH  radiation.  The  corresponding  plots  on  the  right 
side  of  Fig.  9  show  the  determined  averages  of  the  collected  data  points  (solid  lines)  and  the 
sinusoidal  curves  that  were  fitted  to  the  ensemble-averaged  curves  (dashed  lines). 

The  actuator  current,  pressure  and  CH  radiation  are  measured  by  sensors  with  practically  no 
time  delay.  In  contrast,  the  proximity  sensor  that  measures  the  actuator’s  needle  displacement 
introduces  a  significant  phase  lag  and  attenuation  to  the  measured  data.  This  phase  lag  and 
attenuation,  which  are  determined  in  advance  by  comparing  known  input  signals  at  various 
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frequencies  to  the  sensor’s  output,  are  accounted  for  in  the  data  reduction  procedure  by  providing 
the  determined  sinusoidal  signal  with  a  phase  lead  and  a  gain  relative  to  the  calculated  "average" 
curve. 

Figure  9  shows  that  the  measured  CH  radiation  and  pressure  data  form  "clouds”  of  points. 
While  the  CH  radiation  "cloud"  exhibits  a  sine-like  shape,  such  a  pattern  is  not  readily  apparent  in 
the  thick  cloud  of  pressure  data  points.  This  difference  between  the  pressure  and  radiation  data 
is  reflected  in  the  presence  of  high  frequency  components  in  the  ensemble  average  of  the  pressure 
data,  and  the  practically  sinusoidal  shape  of  the  ensemble  average  of  the  CH  radiation  data,  which 
nearly  coincides  with  its  sinusoidal  curve  fit. 

The  above  described  data  reduction  procedure  was  applied  to  the  data  measured  in  the  open 
loop  tests  to  obtain  the  frequency  dependence  of  the  phase  and  gain  of  the  above  described 
variables  and,  in  particular,  the  transfer  function  between  the  needle  displacement,  which  closely 
resembles  the  fuel  injection  modulation,  and  secondary  combustion  process  heat  release 
oscillations. 

The  frequency  dependence  of  the  gain  and  phase  delay  of  the  secondary  combustion  process 
heat  release  oscillations  generated  by  the  fuel  injector  actuator  were  measured  by  the  two,  above 
discussed,  methods  and  the  results  are  compared  in  Figs.  10-a,b.  The  excellent  agreement 
between  the  two  sets  of  data  that  were  determined  using  different  approaches  in  two  different 
combustors  strongly  suggests  that  the  frequency  response  of  the  heat  addition  oscillations 
generated  by  the  developed  fuel  injector  actuator  is  practically  independent  of  the  acoustic 
properties  of  the  combustor. 

Figures  ll-a,b  show  the  measured  frequency  dependence  of  the  phase  and  gain  of  the 
secondary  combustion  process  heat  release  of  two  different  injector/actuator  configurations  that 
consisted  of  each  of  the  developed  injectors,  see  Figs.  7-a,b,  and  the  larger  fuel  injector  actuator 
(i.e.,  configurations  INJ1/ACT2  and  INJ1/ACT2).  The  frequency  dependence  of  the  SCPR  of  the 
INJ1/ACT2  configuration  was  determined  for  various  magnitudes  of  the  actuator’s  needle 
displacement.  Figure  1 1-a  shows  that  all  the  measured  phase  data  can  be  closely  correlated  by  a 
single  curve  that  describes  a  nearly  linear  dependence  of  the  phase  upon  the  frequency.  This  nearly 
linear  frequency  dependence  of  the  phase  suggests  the  presence  of  a  pure  time  delay  between  the 
fuel  injection  rate  modulation  and  the  corresponding  secondary  combustion  process  heat  release 
oscillations.  This  is  a  surprising  result  as  it  was  not  expected  that  an  identical,  pure,  time  delay 
will  be  produced  by  different  injectors  and  excitation  levels. 
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Figure  1 1-b  shows  the  frequency  dependence  of  the  ratio  of  the  amplitudes  of  the  combustion 
process  heat  release  the  actuator’s  needle  displacement  oscillations.  Since  these  results  apparently 
depend  upon  the  amplitude  of  the  needle's  displacement,  this  amplitude  is  indicated  next  to  some  of 
the  data  points  to  provide  means  for  determining  the  dependence  of  the  data  upon  the  needle's 
amplitude.  Figure  1 1-b  indicates  that  the  INJ2  injector  has  a  larger  gain  at  frequencies  below  300 
Hz.  and  smaller  gain  at  higher  frequencies.  It  also  shows  that  the  gain  depends  upon  the  needle's 
displacement  and  is  larger  for  smaller  displacements.  This  result  suggests  that  the  SCPR  is 
controlled  by  nonlinear  processes  whose  potential  origin  is  discussed  below. 

The  developed  fuel  injector  actuator  provides  the  combustor  with  a  fuel  flow  consisting  of  a 
mean  and  an  oscillatory  component.  To  maintain  the  mean  flow  rate  at  a  set  value,  it  is  measured 
and  compared  with  the  set  value.  If  the  two  are  not  equal,  the  mean  flow  rate  controller  sends  a 
slowly  varying  signal  to  the  actuator  that  changes  the  mean  position  of  the  needle  and,  thus,  the 
mean  area  of  the  actuator's  orifice,  in  an  effort  to  provide  the  desired  mean  fuel  flow  rate.  In 
addition  to  the  slowly  varying  signal  from  the  mean  flow  controller,  the  actuator  receives  a  high 
frequency  signal  that  sets  the  needle  into  high  frequency,  back-and-forth,  axial  oscillations  that 
periodically  vary  the  actuator's  orifice  cross  sectional  area,  resulting  in  a  periodic  secondary  fuel 
injection  rate.  As  long  as  the  periodic  variation  of  the  orifice's  area  is  small  and  does  not  change  its 
mean  area,  the  periodic  change  in  the  orifice's  cross  sectional  area  does  not  interfere  with  the  mean 
flow  rate  through  the  actuator.  As  the  amplitude  of  the  actuator's  needle  displacement  oscillations 
increases,  a  threshold  amplitude  that  forces  the  actuator's  orifice  to  momentarily  close  (at  the 
instant  of  maximum"  upward"  needle  displacement,  see  Fig.  6)  is  reached.  As  the  amplitude  of  the 
needle's  displacement  oscillations  increases  beyond  this  threshold  value,  the  mean  values  of  the 
actuator's  needle  displacement  and  mean  flow  rate  increase,  see  Fig.  12.  This  forces  the  actuator's 
mean  flow  rate  controller  to  further  close  the  actuator  orifice,  resulting  in  additional  "truncation"  of 
the  needle's  oscillations.  Since  it  can  be  shown  that  the  amplitude  of  the  actuator’s  orifice  area 
variation  oscillations  cannot  be  larger  than  the  magnitude  the  orifice's  mean  area,  the  amplitude  of 
the  fuel  flow  rate  modulation  cannot  exceed  that  of  the  mean  fuel  flow  rate. 

Figure  13  shows  the  dependence  of  the  ratio  of  the  amplitude  of  the  combustion  process  CH 
radiation  oscillations  and  the  mean  CH  radiation  on  the  needle's  displacement  amplitude  at  a 
frequency  of  540  Hz.  for  three  different  injector/actuator  configurations.  The  behavior  exhibited 
by  the  three  plots  is  in  agreement  with  the  above  discussion.  They  show  that  the  amplitude  of  the 
CH  radiation  oscillations  reaches  a  limiting  value  at  a  given  needle  displacement  amplitude  and 
cannot  increase  beyond  this  limiting  value  as  the  needle  displacement  is  further  increased. 
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To  further  understand  the  results  of  Fig.  13,  we  should  consider  the  two  main  processes  that 
control  the  interaction  between  the  secondary  combustion  process  heat  release  and  needle 
displacement  oscillations.  First,  consider  the  variation  in  the  cross  sectional  area  of  the  actuator's 
orifice,  see  Fig.  6,  in  response  to  the  needle's  displacement  oscillations.  As  discussed  above,  the 
amplitude  of  the  cross  sectional  area  oscillations  cannot  exceed  the  mean  magnitude  of  the  cross 
sectional  area  of  the  orifice.  Consequently,  to  increase  the  amplitude  of  the  fuel  flow  rate 
modulation  would  require  increasing  the  magnitude  of  its  mean  flow  rate.  The  difference  between 
the  behavior  of  the  two  actuators  that  were  used  in  this  study  is  qualitative  described  in  Fig.  14 
where  the  expected  dependence  of  the  amplitude  of  the  secondary  fuel  injection  rate  modulation 
upon  the  amplitude  of  the  actuator's  needle  displacement  is  plotted  for  the  large  and  small  actuators 
operating  with  two  different  mean  flow  rates.  Figure  14  shows  that  since  the  orifice  of  the  larger 
ACT2  actuator  has  a  larger  cross  sectional  area  than  that  of  the  ACT1,  the  ACT2  actuator  provides 
larger  amplitude  fuel  flow  rate  modulations  for  a  given  amplitude  of  the  needle's  displacement 
oscillations  than  the  ACT1  actuator.  It  also  shows  that  both  actuators  can  provide  the  same 
maximum  fuel  flow  rate  oscillations  amplitudes  as  long  as  they  can  provide  the  same  fuel  mean 
flow  rate. 

Another  factor  affecting  the  amplitude  of  the  combustion  process  heat  release  oscillations  is  the 
"efficiency"  of  converting  the  secondary  fuel  injection  rate  modulations  into  combustion  process 
heat  release  oscillations.  Intuitively,  one  would  expect  that  the  ratio  of  the  amplitudes  of  the 
combustion  process  heat  release  and  secondary  fuel  injection  rate  oscillations  would  depend  upon 
the  injector  design.  Examining  the  configurations  of  the  two  investigated  injectors,  see  Figs.  7- 
a,b,  shows  that  injector  ENJ1  supplies  the  secondary  fuel  stream  directly  into  the  "center"  of  the 
primary  combustion  zone  whereas  injector  INJ2  apparently  injects  the  secondary  fuel  upstream  of 
the  primary  combustion  region,  possibly  enabling  the  secondary  fuel  to  fully  or  partially  mix  with 
the  premixed,  primary,  reactants  before  reaching  the  combustion  zone.  This  pre-mixing  process 
apparently  tends  to  decrease  the  magnitude  of  the  attained  secondary  combustion  heat  release 
oscillations.  Consequently,  injector  INJ1  is  expected  to  provide  a  more  "efficient"  conversion  of 
the  fuel  injection  rate  modulation  into  combustion  process  heat  release  oscillations  than  injector 
INJ2. 

On  the  basis  of  the  above  discussion,  one  would  expect  that  configurations  INJ1/ACT2  and 
INJ2/ACT1  will  produce  the  largest  and  smallest  secondary  combustion  process  heat  release 
oscillations  for  a  given  needle  displacement,  respectively,  while  configurations  INJ1/ACT1  and 
INJ2/ACT2  will  produce  heat  release  oscillations  that  fall  between  these  limits.  These  expectations 
are  supported  by  the  results  in  Fig.  13,  where  it  is  shown  that  for  a  given  amplitude  of  the 
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actuators  needle  displacement  the  amplitude  of  the  secondary  heat  release  oscillations  decreased  in 
the  following  order:  INJ1/ACT2,  INJ1/ACT1  and  INJ2/ACT2. 

The  above  discussion  and  Figs.  13  and  14  indicate  that  the  smaller  ACT1  actuator  could 
provide  the  same  maximum  fuel  flow  rate  and  heat  release  oscillations  as  the  larger  ACT2  actuator 
if  it  could  be  provided  with  a  larger  maximum  needle  displacement  amplitude.  This  explains  why 
the  two  configurations  INJ1/ACT1  and  INJ1/ACT2  that  used  the  same  injector  but  different 
actuators  attained  the  same  maximum  secondary  heat  release  magnitudes,  see  Fig.  13.  It  is  also 
interesting  to  note  in  Fig.  13  that  the  maximum  amplitude  of  heat  release  oscillations  that  can  be 
excited  by  the  INJ2  injector  is  significantly  smaller  than  those  excited  with  the  INJ1  injector  even 
though  INJ2  used  the  larger  ACT2  actuator. 
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Closed  Loop  Active  Control  of  Instabilities 


As  stated  earlier,  the  investigated  ACS  is  guided  by  Rayleigh’s  criterion  and  its  objective  is  to 
generate  heat  addition  oscillations  within  an  unstable  combustor  180  degrees  out  of  phase  with  the 
unstable  pressure  oscillations.  This  mode  of  control  will  be  discussed  assuming  that  the  instability 

can  be  damped  by  controlling  its  dominant  mode  only.  In  this  case,  the  observer  only  identifies  the 
dominant  mode  and  the  control  signal  to  the  actuator,  yc,  is  given  by: 

yc  =  G(co)  •  (S  •  sin(ca/  +  4>(co))  +  C  •  cos(atf  +  <p(co)))  (10) 

where  S,  C  and  co  are  the  observed  amplitudes  and  frequency  of  the  dominant  mode  while  G(co) 
and  <J)(co)  are  the  gain  and  phase  that  are  added  by  the  controller  using  data  (in  the  form  of  plots  or 
tables)  obtained  in  the  above  discussed  open  loop  experiments.  Since  Fig.  10-b  shows  that  the 
gain  of  the  heat  addition  oscillations  decreases  with  increasing  frequency,  the  magnitude  of  G(w) 
should  increase  with  frequency  to  assure  that  the  ACS  generates  heat  addition  oscillations  of 
adequate  magnitude  at  higher  frequencies. 


Closed  loop  control  of  the  rocket  motor  setup  shown  in  Fig.  6  was  investigated.  Combustor 
pressure,  control  signal  and  observed  frequency  measured  in  a  typical  test  before  and  after 
activation  of  the  ACS  are  described  in  Fig.  15.  In  this  experiment,  the  development  of  a  large 
amplitude  instability  was  "monitored"  by  the  observer  before  activating  the  control  system  at  t=0.1 
seconds.  The  ACS  was  "conditionally"  activated  by  requiring  the  gain  G(oo)  in  Eq.  10  to  satisfy 
the  following  conditions 


G(co)  =  \ 


0  when  \coobserved  -  o0|>  dx  =  100ft. 


(11) 


[  Gtabk  when  | C0observed  -co0\<d2=  50 Hz. 


where  (Dq  is  the  observed  frequency  of  the  instability  before  activating  the  controller  (see  Fig.  15) 
while  di  and  d2  are  specified  parameters.  These  conditions  require  the  controller  to  "remember" 
the  frequency  of  the  dominant  unstable  mode  before  activating  the  ACS  and  to  respond  only  when 
this  mode  becomes  dominant,  thus  leaving  other  modes  uncontrolled.  This  mode  of  operation 
demonstrates  the  high  flexibility  of  the  developed  control  approach. 
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Figure  15-a  shows  that  this  ACS  practically  damped  the  instability  in  40  milliseconds, 
indicating  that  it  could  damp  rocket  instabilities  before  they  could  seriously  damage  the  engine 
and/or  result  in  mission  failure.  FFT  analysis  of  the  combustor  pressure  oscillations  before  and 
after  the  activation  of  the  ACS  showed  that  the  amplitude  of  the  dominating  mode  was  reduced  by 
26  dB„  which  exceeds  the  performance  of  ACS  investigated  elsewhere  (see  survey  of  results  in 
Table  1  in  Ref.  5).  Figure  15-b  describes  the  time  dependence  of  the  control  signal  to  the  actuator, 
which  is  proportional  to  the  magnitude  of  G(o))  (see  Eqs.  10  and  11)  and,  thus,  depends  upon  the 
observed  frequency,  see  Fig.  15-c.  Figures  15-b,c  show  that  the  control  signal  goes  to  zero  when 
the  370  Hz.  dominant  mode  is  damped  by  the  ACS.  When  this  occurs,  the  observer  identifies  the 
740  Hz.  harmonic  as  the  dominant  mode  and  stops  controlling  the  370  Hz.  mode,  see  Fig.  15-b. 
Consequently,  the  unstable  370  Hz.  mode  starts  growing  again.  When  its  amplitude  exceeds  that 
of  other  unstable  combustor  modes  (e.g.,  the  740  Hz.  mode),  the  observer  identifies  this  mode 
again  as  the  dominant  mode  and  turns  on  the  actuator.  As  shown  in  Figs.  15-b,c,  this  sequence  of 
events  continuously  repeats  itself,  resulting  in  effective  control  of  the  dominant  370  Hz.  mode, 
while  not  destabilzing  other  modes  of  the  combustor. 

Additional  closed  loop  experiments,  not  reported  herein,  have  demonstrated  that  the 
investigated  ACS  can  effectively  control  multi-mode  instabilities. 

Summary  and  Conclusions 

In  closing,  this  report  presents  a  novel  active  control  approach  for  damping  detrimental  rocket 
instabilities.  The  main  contributions  of  this  study  are: 

1 .  the  development  of  an  active  control  approach  based  upon  Rayleigh's  criterion, 

2.  the  development  and  demonstration  of  an  observer  that  can  determine  the  amplitudes, 
phases  and  frequencies  of  a  prespecified  number  of  unstable,  combustor  modes  in  real  time, 

3.  the  development  of  theoretical  models  for  investigating  active  control  of  linear  and 
nonlinear  instabilities  in  rocket  motors, 

4.  the  development  of  a  fuel  injector  actuator  that  utilizes  a  magnetostrictive  material  to 
modulate  a  secondary  fuel  injection  rate  over  a  0-1,200  Hz.  range,  which  is  wider  than  that  of  any 
known  injector, 

5 .  the  demonstration  that  the  developed  fuel  injector  actuator  can  excite  significant  reaction  rate 
heat  release  oscillations  within  the  combustor. 
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6.  the  development  of  two  different  approaches  for  determining  the  frequency  dependence  of 
the  response  of  the  secondary  combustion  process  generated  by  the  fuel  injector  actuator,  and 

7 .  the  demonstration  that  closed  loop  application  of  the  investigated  ACS  significantly  (e.g.  by 

26  dB.)  and  rapidly  (e.g.,  within  40  milliseconds)  damps  a  rocket  motor  instability  without 
destabilization  of  any  stable  modes. 

It  is  believed  that  these  findings  provide  a  foundation  for  guiding  the  development  of  ACS  for 
unstable  rockets  and  other  combustors. 
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Figure  1.  A  schematic  of  the  developed  ACS. 
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Figure  2-c.  Observed  time  dependence  of  the  secondary  mode  of  the  instability. 


Figure  2-d. 


0.3 


■Measured 

Reconstructed 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


Reaction  rate  model  (Global  Kinetics) 

368] 

®/  =  ®ch.  =  -2.15xlOl6[pc„J  [paj  r  j 

Heat  release 
q^-QcaiU, 

Weighting  function 

x>xf 

=  0  x<xf 

Umb  Weighted  average  of  the  state  vector 

s},s2,s3  Source  terms  due  to  secondary  mass,  momentum 
-and  heat  addition 

Sj  Source  term  due  to  secondary  (controller)  fuel  addition 

q  Heat  release  due  to  chemical  reaction  (j3  contains  q) 

Qai  Calorific  value  of  the  fuel  (methane) 

xf  Flame  holder  location 

lmiI  Characteristic  mixing  length 

Figure  3.  Conservation  equations  used  to  model  the  behavior  of  an  unstable  gaseous  rocket 
motor. 
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Figure  4.  Prediction  of  the  combustor  response  to  open  loop  excitation  by  oscillatory  fuel 
injection. 
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Figure  5.  Results  of  numerical  simulation  of  an  unstable  combustor  response  to 
closed  loop  active  control  by  means  of  an  oscillatory  secondary  fuel 
injection  using  real  time  observation  of  the  dominant  mode. 
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Time  dependence  and  spectra  of  the  combustor  pressure  and  combustion 
region  CH  radicals  radiation  measured  in  an  open  loop  experiment. 
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Figure  9. 


Measured  and  ensemble  averaged  data  describing  the  actuator  current, 
needle  displacement,  combustor  pressure  and  CH  radicals  radiation  obtained 
in  the  open  loop  experiment  whose  results  are  presented  in  Fig.  8. 
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Phase  difference  between  heat  release  and  actuator  displacement  oscillations 
for  two  different  injector/actuator  configurations. 
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Figure  12.  Qualitative  description  of  the  effect  of  mean  flow  control  upon  fuel  flow  modulation 
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The  dependence  of  the  ratio  of  oscillating  to  mean  heat  release  radiation 
upon  the  displacement  of  the  actuator's  needle  amplitude. 


Figure  14-a.  Qualitative  description  of  the  dependence  of  the  fuel  modulation  amplitude 
upon  the  needle  displacement  amplitude  for  two  different  actuators  and 
mean  flows. 


Figure  14-b.  Qualitative  description  of  the  dependence  of  the  heat  release  amplitude  upon 
the  amplitude  of  the  fuel  flow  modulation  for  two  different  injectors. 
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Appendix  A 

Extension  of  the  Theory  of  the  Observer  for  Simultaneous 
Identification  of  an  Arbitrary  Number  of  Modes  in  the  Input  Signal 

Prepared  by 

Nikos  Markopoulos  and  Yedidia  Neumeier 


1.  Introduction 

The  frequency  adaptation  formula  used  so  far  in  our  observer  is  based  on  a  single,  sinusoidal  input. 
If  we  know  a  signal  at  each  instant  of  time,  and  if  we  also  know  that  the  signal  is  indeed  a  sinusoid,  then, 
this  frequency  adaptation  formula  supplies  us  with  the  exact  frequency  of  the  sinusoid,  without  the  need 
for  any  iterations.  In  practice,  the  need  for  iterations  arises  because  of  at  least  two  reasons.  First,  the  real 
pressure  signals  that  we  encounter  in  a  combustion  chamber  vary  with  time  in  a  manner  that  is  far  more 
complex  than  that  of  a  single  sinusoid.  Second,  the  frequency  domain  characteristics  of  such  signals  are 
not  stationary,  but  change,  at  best  slowly,  with  time.  We  may  be  able  to  improve  this  situation  at  least 
partially  by  generalizing  our  frequency  adaptation  formula,  so  that  it  is  accurate  for  input  signals  that  are 
composed  of  an  arbitrary  number  of  sinusoids.  Then,  by  accounting  for  a  large  enough  number  of 
sinusoids,  we  hope  to  increase  the  accuracy  and  ease  with  which  we  capture  the  characteristics  of  real 
signals  exhibiting  a  discrete  energy  spectrum.  This  in  turn  may  make  our  observer  more  robust  and  reduce 
some  of  our  worries  about  convergence.  This  generalization  of  the  frequency  adaptation  formula  is  given 
here  in  two  steps.  First,  in  Section  2,  the  formula  is  generalized  to  an  input  signal  that  is  equal  to  a  sum  of 
two  sinusoids.  Then,  with  the  insights  gained  in  Section  2,  the  formula  is  generalized  in  Section  3  to  an 
input  signal  that  is  equal  to  a  sum  of  an  arbitrary  number  of  sinusoids. 

2.  Extending  the  Theory  of  the  Observer  to  two  modes 

2.1  The  input  signal 

Consider  a  signal  f(t)  given  by: 

f(t)  =  f1(t)  +  f2(t)  (1) 

where  fj(t)  and  f2(t)  are  two  sinusoids: 

f  l(t)  =  Cisin£21t  +  D1cosQ1t  (2) 

f2(t)  =  C2sinfl2t  +  D2cosQ2t  (3) 

Assume  that: 

(i)  We  know  f(t)  at  each  instant  of  time  t. 

(ii)  We  know  that  f(t)  is  a  sum  of  two  sinusoids. 

The  problem  that  we  will  try  to  solve  is  the  following:  Based  on  the  information  given  in  (i),  (ii), 
determine  the  characteristics  of  f(t),  that  is,  the  frequencies  Qj  and  ft2,  and  the  amplitudes  Clt  Dj,  C2, 
and  D2. 

2.2  The  observer 

To  solve  the  problem,  in  analogy  with  the  one-mode  case,  we  set  up  an  “observer”  signal  F(t)  in  the 
following  way:  F(t)  is  equal  to  the  sum: 
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with  Fj(t)  and  F2(t)  given  by: 


F(t)  =  F1(t)  +  F2(t) 

(4) 

Fi  (t)  =  A](t)sino)it  +  B!(t)coso)it 

(5) 

F2(0  =  A2(t)sinco2t  +  B2(t)cosco2t 

(6) 

Here  coj  and  w2  are  two  arbitrary  frequencies  that  we  select  in  advance  for  our  observer.  The  time 
functions  Aj(t),  Bj(t),  A2(t),  and  B2(t)  are  given  by: 

t 

^1(0  =  ^-  J[f(T)-F2(T)]sinco1T  dx 


t-T, 
t 


Bi(0  =  ^r-  J[f(x)-F2(x)]cosco1x  dx 


t-T, 

t 


a2(0-^-  J[f(x)-Fj(x)]sin(o2x  dx 


(7) 


(8) 


(9) 


t-T? 


I 

b2(0  =  —  J[f(x)-F1(x)]cosco2x  dx 


t-T, 


(10) 


In  Eqs.  (7-10)  T1=2jt/o)1  and  T2=27t/o)2  are  the  periods  corresponding  to  the  observer  frequencies.  Using 
some  well-known  trigonometric  identities,  we  can  completely  eliminate  Aj(t),  Bj(t),  A2(t),  and  B2(t)  and 
work  directly  with  Fj(t)  and  F2(t).  Combining  Eqs.  (5),  (7),  and  (8)  one  obtains: 


l 

Fi(0  =  J  [f(x)  -  F2(x)]cost0!(t  -  x)  dx 


(11) 


t-T, 


Similarly,  combining  Eqs.  (6),  (9),  and  (10)  results  in: 


l 

f2(0  =  -^-  J[f(x)-F1(x)]cosco2(t-x)  dx 


(12) 


t-T, 
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Our  observer  is  completely  described  by  Eqs.  (4),  (1 1),  and  (12).  If  one  thinks  of  F^t)  and  F2(t)  as 
the  estimates  of  fj(t)  and  f2(t)  respectively,  then,  from  Eq.  (11),  the  input  for  evaluating  the  estimate  of 
fl(t)  is  f(t)  minus  the  estimate  of  f2(t).  Likewise,  from  Eq.  (12),  the  input  for  evaluating  the  estimate  of 
f2(t)  is  f(t)  minus  the  estimate  of  fj(t).  The  symmetry  is  complete. 

For  the  particular  f(t)  given  by  Eqs.  (1-3),  the  solution  of  Eqs.  (1 1-12)  can  be  expressed  as: 


Fi(t)=Fn(t)  +  F12(t) 

(13) 

where  Fn(t)  and  Fj2(t)  are  two  sinusoids: 

Fj  i  (t)  =  Kj  i  sin  Qjt  +  Ljj  cosftjt 

(14) 

Fj2(t)  =  Kj2  sinQ2t  +  Lj2  cosQ2t 

(15) 

and: 

F2(t)  -  F2j(t)  +  F22(t) 

(16) 

where  F2j(t)  and  F22(t)  are  the  two  sinusoids: 

F21(t)  =  K21  sinQjt  +  L2j  cosQjt 

(17) 

F22  (0  =  K22  sin  Q2t  +  L22  cos  Q2l 

(18) 

The  amplitudes  Ky  and  Ly  are  independent  of  t  They  depend  in  quite  a  complex  manner  on  the 
frequencies  Oj  and  Q2  and  amplitudes  Cj,  Dj,  C2,  and  D2  of  the  input  signal  f(t),  plus  the  frequencies  coj 
and  co2  of  our  observer.  To  Find  this  dependence  one  can  substitute  from  Eqs.  (13-18)  into  Eqs.  (11-12) 
and  determine  Ky,  Ly  by  equating  the  coefficients  of  terms  in  sinQjt,  cosQjt,  sinQ2t,  and  cosQ2t.  This 
step  is  very  tedious  and  time  consuming  and  won’t  be  given  here  explicitly. 

2.3  Perfect  reconstruction  of  the  input 

When  our  observer  operates  at  the  “correct”  frequencies,  that  is,  when  and  co2=Q2,  then,  one 

can  show  that  the  cross-coupling  sinusoids  F12(t)  and  F21(t)  are  identically  zero,  while  F1(t)=F11(t)=f1(t), 
and  F2(t)=F22(t)=f2(t).  In  this  case  fj(t)  is  perfectly  reproduced  by  Fi(t),  f2(t)  is  perfectly  reproduced  by 
F2(t),  and  the  input  signal  f(t)  is  perfectly  reproduced  by  the  observer  signal  F(t).  Thus,  in  order  to  find 
the  amplitudes,  phase  shifts,  and  frequencies  of  the  input  signal  f(t),  all  we  have  to  do  is  determine  first,  in 
some  way,  the  frequencies  Qi  and  Q2.  After  that,  we  can  simply  update  the  frequencies  of  our  observer  to 
and  £22.  Then,  the  signals  Fj(t)  and  F2(t)  produced  by  our  observer  will  be  the  same  as  the  sinusoids 
fl(t)  and  f2(t)  respectively. 

2.4  Determining  the  correct  frequencies 

The  reason  that  the  determination  of  and  Q2  is  not  straightforward  is  that  our  observer  supplies 
us  only  with  the  signals  Fjlt)  and  F2(t).  That  is,  we  have  no  knowledge  of  the  components  Fn(t),  F]2(t), 
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F2i(t),  and  F22(t)  individually.  Similarly,  we  have  no  knowledge  of  the  components  fj(t)  and  f2(t)  of  f(t) 
individually.  All  we  know  is  their  sum  f(t).  Keeping  this  in  mind,  the  frequencies  ft]  and  ft2  can  be  found 
by  extending  the  method  applied  previously  to  a  one-mode  signal  f(t)  in  the  following  way:  First,  one  can 
examine  the  explicit  expressions  for  the  amplitudes  Kij  and  Lij  (not  given  here)  and  show  that  the  time 


derivatives  of  the  functions  Fjj(t)  can  be  written  as: 

2 

Fll(0  —  /“2  2\  [f  1  (t) - F21  ( t) - f i(t - T, ) + F21(t - Tj )]  (19) 

7t^ft]  -  CD]  J 
2 

f12(t)  =  -7X22\  [f 2 (0 - F22 (t ) - f 2 (t - T, )  +  F22 (t - T] )]  (20) 

7C^ft2-(0]  J 

2 

*21(0 - rai\  [fl(')-Fll(l)-fl(t-T2)  +  Fn(t-T2)]  (21) 

2 

F22(0-  ,“2,  2  2,  [f2(0 - F12(t) - f2(t - T2)+  Fl2 (t - T2)]  (22) 

7t^ft2“  C02  J 

Then,  combining  Eqs.  (19-20),  and  using  Eqs.  (1),  (16)  one  obtains: 


Fj  J  (t) + f  F12  (t)  =  [f  (t>  -  F2  (‘ )  -  f  (t  -  T! )  +  F2  (t  -  Tt)]  (23) 

£2i  ;  “2  ;  n 

Similarly,  combining  Eqs.  (21-22),  and  using  Eqs.  (1),  (13)  results  in: 

22«)  =  ^  [f(0  -  F,(t)  -  f  (t  -T2)  +  F,(t  -  T2)]  (24) 

Denoting  the  terms  in  the  brackets  above  by  M](t)  and  M2(t): 

M](t)  =  f  (t)  -  F2(t)  -  f(t  -  T] )  +  F2(t  -  T])  (25) 

M2(t)  =  f(t)-F](t)-f(t-T2)  +  F](t-T2)  (26) 

and  using  Eqs.  (13),  (16),  one  can  write  Eqs.  (23),  (24)  as: 
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(27) 


2  2 

T^-FiiO)  +  7^12(0  =  F,(t)  -  ^  Ml(t) 

£2j  122  71 

2  2 

^■F2l(0  +  ^-F22(0  =  F2(t)-2iM2(t)  (28) 

12j  122  71 

Defining  two  auxiliary  variables  by: 


Ni(t)  =  F^t)  -  Mj(t) ,  N2(t)  =  F2(t)  -  ^  M2(t)  (29) 

7t  It 


and  using  again  Eqs.  (13),  (16),  one  can  further  manipulate  Eqs.  (27)  and  (28)  into  the  form: 


nF*l(‘)+ 


1  1 


w 


'12  (0  «  M 

(Of 


(30) 


f 


w& 


2  ) 


F2l(t)  +  7^F2(t)  = 

“2 


N2(l) 

2 

<02 


(31) 


We  want  to  solve  Eqs.  (30-31)  for  £2j,  fi2’  but  the  problem  is  that  we  don’t  know  the  cross¬ 
components  Fj2(t),  F2i(t)  (see  Eqs.  (13-18)).  If  the  observer  frequencies  (Oj,  o>2  are  close  enough  to  the 
input  frequencies  respectively,  then,  by  examining  the  corresponding  expressions  for  the 

amplitudes  Ky,  Ly,  of  Fy(t)  (not  given  here  explicitly),  one  can  show  that  K12,  L12,  and  K2\,  L2i,  are 
much  smaller  compared  to  Kjj,  Ln,  and  K22,  L22>  respectively.  In  such  a  case,  it  is  safe  to  neglect  the 
time  derivatives  of  Fj2(t),  F21O)  in  Eqs.  (30-31),  and  solve  these  equations  approximately  for  Qj,  Q2  as: 


q2=  «©fFl(t) 

1  n  Fj  (t)  -  coj  [f (t )  -  F2  ( t) -  f (t  -  Tj )  +  F2  ( t  -  Tj ) 


(32) 


q2  _  7t  <02F2(t) 

7t  F2(t)  -  (02  [f (t)  -  Fj  (t)  -  f  (t  -  T2 )  +  Fj  (t  -  T2 )] 


(33) 


2.5  Comparison  with  the  old  frequency  adaptation  scheme 

Before  proceeding  further,  it  is  interesting  to  see  the  connection  between  Eqs.  (32-33)  and  the 
frequency  adaptation  used  in  our  present  day  observer.  To  make  this  connection,  consider  first  an  input 
signal  f(t)  which  is  a  single  sinusoid  with  frequency  £2.  If  F(t)  is  the  corresponding  (single)  observer 
signal  that  we  use  to  reproduce  f(t),  and  if  ©  is  the  frequency  of  our  observer,  then  £2  can  be  found  from 
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the  old  adaptation  formula,  derived  by  Dr.  Neumeier,  who  worked  with  the  quantities  A(t)  and  B(t)  rather 
than  F(t).  Using  our  notation  this  old  adaptation  formula  can  be  written  as: 


r_2_  7tto2F(t) 

7t  F(t)  —  co  [f (t)  —  f (t  —  T)] 


(34) 


Eq.  (34)  supplies  us  with  the  exact  value  of  Q,  as  long  as  f(t)  is  a  single  sinusoid.  If  we  now  consider  an 
input  signal  f(t)  which  is  a  sum  of  n  sinusoids,  the  intuitive  way  to  generalize  Eq.  (34)  would  be: 


Fj(t) 


7C  Fj(t)-COi  > 

"  *  1  V  V 

n 

f(t)  —  f^-TiJ-^FjW-Fjtt-Ti)]. 

j=w 

i  =  1,2, 


n 


(35) 


Clearly,  Eqs.  (32-33)  correspond  to  the  special  case  n=2  in  Eq.  (35).  What  is  done  in  our  present  day 
observer  corresponds  basically  to  using  Eq.  (35)  in  an  iterative  scheme  to  determine  the  frequencies  Qj  of 
an  input  signal.  The  iterations  are  needed  partly  because  Eq.  (35)  is  only  approximate,  since  it  ignores  the 
coupling  terms  Fjj(t),  and  partly  because  real  signals  may  have  characteristics  that  change  with  time,  like 
component  frequencies,  amplitudes,  phase  shifts,  etc. 


2.6  Accurate  frequency  adaptation  for  two  modes 

We  can  try  to  improve  this  situation  by  obtaining  a  more  accurate  way  of  determining  the  component 
frequencies  Qj  for  input  signals  composed  of  n  sinusoids.  We  will  first  do  this  for  the  case  n=2,  then 
generalize  it  to  an  arbitrary  n. 

For  n=2  we  have  to  reconsider  the  exact  Eqs.  (30-31).  These  equations  are  valid  at  any  instant  of 
time  t,  in  particular,  they  are  valid  for  the  time  instants  t-At,  and  t+At,  where  At  can  be  taken  as  a  known, 
strictly  positive,  and  apart  from  that  arbitrary  time  increment.  Explicitly,  at  t-At  one  has: 


— TF1(t-At)  + 
i2l 


of  of 


Fi2(t-At) 


Nj(t- At) 
- 3 - 


(36) 


1 


1 


Qf  nj 


F2l(c-dt)  +  -4.F2(t-At)=  -N2<t  i  Al) 


(37) 


and,  at  t+At  one  similarly  has: 


A-7 


jp-F^t  +  AO  +  f^-jp-  jF12(t  +  At)  = 

i2i  ^£22  £2j  )  COj 


%~k  M‘+*>+i***<*+*>-  ^ 


From  Eqs.  (15),  one  can  now  show  that: 


F12(t  -  At)  =  F12(t)  cos(Q2At)  +  n2  F12(t)  sin(Q2At) 
Fi2(t  + At)  =  F12(t)  cos(f22At)  -  £22F12(t)  sin(£22At) 


Similarly,  from  Eqs.  (17),  one  can  show  that: 


F2j  (t  -  At)  =  F2i  (t)  cos(f2jAt)  +  F2j  (t)  sin(QjAt) 

F2](t  + At)  =  F2j (t)  cos(f2jAt)  -  Qj F2j(t)  sin(QiAt) 

Substituting  back  from  Eqs.  (40-43)  into  Eqs.  (36-39)  we  get: 

■~jF1(t-At)+  73-77 j  [^12(0  cos(Q2At)  +  &2F12(t)  sin(£22At)]  =  —■-•7—" 
£2i  ^£22  £2j  J  COj 


3-73-|[f21(t)cos(QiAt)  +  n1F2i(t)sin(n,At)]+4j.t2(t-At)=  ML-^1  (45) 

»2j  i22  J  s22  C02 


-i-F1(t  +  At)+  73-73-  |[F12(t)cos(f22At)-£22F12(t)sin(£22At)]  =  (46) 

i2l  v*22  £2!  J  COi 


73-7^  P21W  cos(QiAt)  -  O,  F21(t)  sin(Q1At)]  +  -4-F2(t  + At)  =  (47) 

£2f  )L  coj 


Equations  (44-47),  together  with  Eqs.  (30-31)  constitute  a  set  of  six  equations  in  six  unknowns.  The 
unknowns  are  the  four  quantities  Qj,  Q2,  F12(t),  F21(t),  and  the  time  derivatives  of  F12(t)  and  F21(t). 
From  Eqs.  (30),  (44),  and  (46),  eliminating  F12(t)  and  its  time  derivative  one  obtains: 
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r>2_m2[  2F1(t)cos(fl2At)-F1(t  +  At)-F1(t-At)1 

1  1  [2  N1(t)cos(D2At)  -  Ni(t  +  At)  -  Nj(t  -  At)  J  (48) 

Similarly,  from  Eqs.  (31),  (45),  and  (47),  eliminating  F2j(t)  and  its  time  derivative  results  in: 

q2  =  J  S^Wcos^Atj-Fza  +  AO-^a-At)  1 

2  2  _  2  N2  (t)  cos(Qj At)  -  N2  (t  +  At)  -  N2  (t  -  At)  ' 

Equations  (48-49)  must  now  be  solved  for  the  two  unknown  frequencies  and  Q2.  In  general,  this 
requires  the  application  of  a  numerical  iterative  scheme.  If  we  select  a  very  small  time  increment  At,  say, 
equal  to  the  integration  time  step  for  the  observer  signals  Fj(t)  and  F2(t),  then,  we  can  presumably 
approximate  the  cosines  in  Eqs.  (48-49)  as: 


cos(DjAt)  =  1 - i 


n?At2 


cos(f22At)  =  1 - - 


Qo  At2 


In  this  case  Eqs.  (48-49)  can  be  solved  analytically.  If  At  is  small  enough  so  that  the  cosines  can  be 
approximated  by  unity,  then,  Eqs.  (48-49)  decouple  and  supply  us  with  the  solution: 

r>2  2|~  2Fi(t)~Fl(t  +  At)-Fi(t-At)  1  ... 

^1  o  xt  /+\  xt  ,  a  xt  /  ax’  as  At  — >  0  (51) 

_2Ni(t)-Nj(t  +  At)-Nj(t-  At) 


q2  _  m2[  2F2(t)-F2(t  +  At)-F?(t-At)  ' 
2|_2N2(t)-N2(t  +  At)-N2(t-At)  ’ 


as  At  — »  0 


for  the  frequencies  and  Q2.  This  solution  tends  to  the  exact  solution  as  At  tends  to  zero.  In  practice,  for 
a  finite,  small  At  the  solution  expressed  in  Eqs.  (51-52)  is  either  close  enough  to  the  exact  solution,  or  if 
not,  it  can  be  used  as  a  first  guess  for  the  iterative  solution  (at  a  fixed  t)  of  Eqs.  (48-49). 

3.  Extending  the  Theory  of  the  Observer  to  n  modes 

3.1  The  input  signal 

With  the  two-mode  case  in  the  background,  the  general  case  of  n  modes  now  becomes  more 
transparent.  Consider  again  a  signal  f(t)  given  by: 


n 

f(0  =  2/i(t) 


where  fj(t)  are  individual  sinusoids: 
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i  =  1,2,  ....,n 


(54) 


fj(t)  =  CiSin£2jt  +  Dj  cosQjt , 

Assume  again  that: 

(i)  We  know  f(t)  at  each  instant  of  time  t. 

(ii)  We  know  that  f(t)  is  a  sum  of  n  sinusoids,  where,  n  is  known. 

Again,  based  on  the  information  given  in  (i),  (ii),  we  would  like  to  determine  the  frequencies  Qj,  and 
the  amplitudes  Cj,  Dj  of  f(t). 


3.2  The  observer 

To  solve  the  problem,  we  set  up  again  an  observer  signal  F(t)  which  is  equal  to  the  sum: 

n 

FCO-XW 

i=l 

where,  in  analogy  with  Eqs.  (11-12),  F,(t)  are  given  by: 


coso)i(t-x)  dz , 


i  =  1,2, ....,n 


(55) 


(56) 


The  amplitudes  Ky  and  Ly  are  independent  of  L  As  in  the  two-mode  case,  they  depend  on  the 
frequencies  and  amplitudes  Cj,  Dj,  of  the  input  signal  f(t),  plus  the  frequencies  cOj  of  the  observer.  To 
guess  this  dependence  one  substitutes  from  Eqs.  (57-58)  into  Eqs.  (56)  and  determines  Ky,  Ly  by 
equating  the  coefficients  of  terms  in  sinf2jt  and  cosQjt. 


3.3  Perfect  reconstruction  of  the  input 

In  analogy  with  the  two-mode  case,  one  can  show  that  an  observer  operating  at  the  frequencies 
rnpQi  reproduces  the  input  signal  perfectly.  In  this  case  all  the  cross-coupling  sinusoids  Fy(t),  i*j,  are 
identically  zero,  while  Fii(t)=fi(t).  Thus,  to  find  the  characteristics  of  the  input  signal,  we  must  again  first 
determine  the  input  frequencies  £2j. 
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3.4  Determining  the  correct  frequencies 

Before  determining  the  frequencies  we  have  to  note  that  we  know  only  f(t)  and  the  signals  Fj(t) 
supplied  to  us  by  our  observer.  We  have  no  knowledge  of  the  individual  components  Fjj(t)  of  the 
observer,  or  the  individual  components  fj(t)  of  the  input.  Luckily,  Eqs.  (19-22)  which  led  to  the 
determination  of  Qj  in  the  two-mode  case  afford  a  straightforward  generalization  for  the  present  case.  By 

examining  the  explicit  expressions  for  the  amplitudes  Kij  and  Lij  (not  given  here)  one  can  show  that  the 
time  derivatives  of  the  functions  Fjj(t)  can  be  written  as: 

2  f  n 

CO;  Qf  rn  _ 

Fii(‘)=»[fif-m?)  fj(l)-fj('  — T‘)  (59) 

1  [  k=l,k*i 

Rearranging  the  above  expression  and  summing  over  the  index  j  one  obtains: 

r 

_ n  ^  /  2  2  A  n  ^  n  n 

j"1  l  j=l  j=l  k=l,k/i 

The  sum  on  the  left-hand -side  can  be  simplified  by  using  Eq.  (57).  The  first  sum  on  the  right-hand-side 
can  be  simplified  by  using  Eq.  (53).  For  the  double  sum  on  the  right-hand-side  one  first  uses  Eq.  (57)  and 
sums  with  respect  to  j,  then  the  dummy  index  k  is  changed  to  j.  After  rearranging  terms,  Eq.  (60)  can  be 
written  in  the  equivalent  form: 


n 

I 


T  Fij(t)  =  Fi(t)-^  f(t)-f(t-Ti)  - 


n 


Just  as  in  Eqs.  (25-26),  defining  the  auxiliary  quantities  Mj(t)  by: 


n 

Mj(0=  f(t)-f(t-Tj)-  yVftJ-Fjd-Ti)], 


j=l,  j/i 


one  can  write  Eq.  (61)  as  (see  Eqs.  (27-28)): 


n 

2 


$  Fjj(t)  =  Fj (t)  — i-Mi(t) 


i  =  1,2, . ,n 


Then,  using  the  definitions  (see  Eqs.  (29)): 


A-ll 


Nj(t)  =  Fj(t)  -  —  Mj(t) ,  i  =  1,2 . ,n 

7t 

and  Eq.  (57),  one  can  write  Eq.  (64)  in  a  form  completely  analogous  to  Eqs.  (30-31): 

n 

eff'(0  +E 


_  Nj(t) 


CO? 


i  =  1,2, . ,n 


(64) 


(65) 


To  solve  the  above  n  equations  for  £2j  we  first  need  to  have  sufficient  information  to  determine  the  cross¬ 
coupling  sinusoids  Fy(t).  Extending  the  procedure  followed  in  section  2.6,  we  can  select  an  integer  m,  and 
write  equations  identical  to  Eqs.  (65),  valid  at  the  2m  instants  of  time  t-kAt,  t+kAt,  with  k=l,2,....,m.  We 
thus  obtain  an  additional  2mn  equations  (just  like  Eqs.  (36-39)  in  the  two-mode  case).  In  each  one  of  these 
equations,  just  as  in  Eqs.  (40-42),  we  use  the  identities: 

Fy  (t-kAt)  =  Fjj(t)  cos(k£2jAtj  +  Qj  Fy(t)  sin(kQjAt)  (66) 

Fij(t  +  kAt)  =  Fjj(t)  cos(k£2jAt)  -  Qj  Fy(t)  sin(kQjAt)  (67) 


Then,  we  have  a  total  of  2n2-n  unknowns.  The  unknowns  are  the  n  frequencies  Qj,  the  n2-n  quantities 
Fy(t),  and  the  n2-n  derivatives  of  Fy(t).  We  also  have  2mn+n  equations,  n  of  them  given  by  Eqs.  (65), 
and  2mn  of  them  given  by  (see  Eqs.  (44-47)): 


n  ,  \ 

■— 5fFi(t  — kAt)  +^)^  I  —  -^2  ^Fy(t)  cos^kf2jAtj  +  fij  Fjj(t)  sin^kQjAtjj  = 


Nj  (t-kAt) 
1 - 


CD 


(68) 


n  ,  \ 

^2  Fi (t  +  kAt)  +^ 1 7^2  ~  ^2  [Fjj(t)  cos(kQjAt)  -  Dj  Fy(t)  sin(kfijAt)]  = 
j=l,j*i  1  1 


Nj  (t  +  kAt) 
CD? 


(69) 


Sufficient  information  for  a  solution  is  thus  provided  when  2mn+n=2n2-n.  This  determines  the  value  of  m 
as  m=n-l.  Adding  Eqs.  (68-69)  one  obtains: 

n 


Q 


77 


Fy(t)  cos^kQjAtj  = 


_  Nj(t  +  kAt)  +  Nj(t-kAt)  Fj(t  +  kAt)  +  Fj(t-kAt) 


2CDj 


2  Qf 


(70) 
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For  m=n-l,  Eqs.  (70)  constitute  a  set  of  mn=n2-n  equations  that  can  be  used  to  eliminate  the  n2-n  time 
derivatives  of  Fjj(t).  The  resulting  expressions  for  these  time  derivatives  can  then  be  substituted  into  Eqs. 
(65),  supplying  us  with  a  system  of  n  equations  in  the  n  frequencies  Qj.  Just  as  in  the  two-mode  case,  a 
complete  solution  for  the  frequencies  Qj  can  be  obtained  when  the  time  increment  At  tends  to  zero.  In  this 
case  all  the  cosines  can  be  safely  replaced  by  unity,  and  Eq.  (70)  becomes: 
n 

j=l,j*i 

Substituting  from  Eq.  (71)  into  Eq.  (65),  and  taking  k=l,  results  in  a  set  of  n  decoupled  equations  which 
can  be  solved  for  the  frequencies  Qj  as: 


_  N^t  +  kAt)*  Nj(t-kAt)  Fj  (t  +  k  At)  -f  F;  (t  -  k  At) 


2  of 


2  Clf 


as  At  0  (71) 


2Fi(t)-Fi(t  +  At)-Fj(t-At) 
2Nj(t)-  Nj(t  +  At)  —  Nj  (t  —  At) 


as  At— >0,  i  =  1,2, . ,n 


(72) 


This  solution  is  qualitatively  identical  to  the  one  obtained  for  the  two-mode  case  (compare  with  Eqs.  (51- 
52)),  and  tends  to  the  exact  solution  for  fij  as  At  tends  to  zero.  In  practice,  for  a  finite,  small  At  the 
solution  expressed  in  Eqs.  (72)  is  either  close  enough  to  the  exact  solution,  or  if  not,  it  can  be  used  as  a 
first  guess  for  the  iterative  solution  (at  a  fixed  t)  of  Eqs.  (65),  (70). 
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Abstract 

This  paper  discusses  the  feasibility  of  high  frequency  nonlinear  vibrational  control. 
Such  control  has  the  advantage  that  it  does  not  require  state  measurement  and 
processing  capabilities  that  are  required  in  conventional  feedback  control.  Bellman, 
Bentsman,  and  Meerkov  (1986)  investigated  nonlinear  systems  controlled  by  lin¬ 
ear  vibrational  controllers  and  proved  that  vibrational  control  is  not  feasible  if  the 
Jacobian  matrix  has  a  positive  trace.  This  paper  extends  previous  work  to  include 
nonlinear  vibrational  controllers.  A  stability  criteria  is  derived  for  nonlinear  systems 
with  nonlinear  controllers,  and  it  is  shown  that  a  nonlinear  vibrational  controller 
can  stabilize  a  system  even  if  the  Jacobian  matrix  has  a  positive  trace. 
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1  Introduction 


This  paper  discusses  the  feasibility  of  applying  open  loop  control  in  the  form  of  high  frequency  vibrational  control 
to  engineering  systems.  Such  control  may  be  applied  in  cases  where  closed  loop  control  is  impractical  and  has  the 
advantage  that  it  does  not  require  costly  sensing  and  computing  capabilities.  Vibrational  control  is  applied  by 
oscillating  an  accessible  system  component  at  low  amplitude  and  high  frequency  (relative  to  the  natural  frequency 
of  the  system).  For  example,  an  inverted  pendulum  can  be  stabilized  by  vertically  oscillating  the  pendulum  pin 
at  a  sufficiently  high  frequency  and  low  amplitude.  Let  us  examine  the  case  of  the  pendulum  in  more  detail.  The 
vertically  oscillated  pendulum  is  described  by  the  following  nonlinear  differential  equation, 

=  *2,  (1.1) 

X2  =  Csin(xi)  -  Bx2  +  aw2Ds\n(xi)s\n(wt),  (1.2) 

where  xi  is  the  angular  displacement  measured  from  the  inverted  equilibrium  point,  X2  is  the  angular  velocity,  B, 
C  and  D  are  positive  physical  constants,  and  a  and  w  are  the  amplitude  and  frequency  of  the  applied  vibration, 
respectively.  In  this  example,  the  control  input  is  the  applied  vibration,  which  is  given  by  asin(utf).  Note  that 
the  amplitude  and  frequency  of  the  control  input  are  constant  and,  therefore,  independent  of  the  state  of  the 
system.  Since  there  is  no  sensing  or  computation  involved,  this  is  a  form  of  open  loop  control.  However,  (1.2) 
involves  a  feedback-like  term  w1  D sin(xi),  which  occurs  naturally  as  a  result  of  the  moment  arm  sin(xi)  between 
the  vertically  oscillating  pendulum  pin  and  the  center  of  mass  of  the  pendulum.  Consequently,  the  feedback 
u,2/?sin(xi)  is  naturally  occurring . 


Vertically 

vibrated 

pendulum 


Horizontally 

vibrated 

pendulum 


Since  the  naturally  occurring  feedback  i^Dsin^)  in  (1.2)  is  of  the  same  form  as  Csinfo),  we  can  view  this 
form  of  control  as  a  variation  of  the  parameter  C;  that  is, 

*1  =  *2,  (1.3) 

X2  =  [C  +  oti>2Z)sin(t^)]sin(ii)  -  Bx7.  (1.4) 

Linearization  of  the  above  system  yields 

0  1  xi 

C  +  aw2Dsin{xvt)  -B  x2 

which  is  of  the  form 

x  =  (.4  +  B(t)]x,  (1.6) 

where  x  is  a  vector,  A  is  a  constant  matrix  and  B(t)  is  a  time- varying  matrix.  In  the  linear  model  (1.6),  vibrational 
control  appears  as  a  variation  of  parameters,  where  the  parameters  of  the  matrix  A  are  varied  by  B(t).  This 
is  the  model  investigated  by  Bellman,  Bentsman  and  Meerkov  [1],  However,  there  is  no  reason  to  assume  that 
vibrational  control  can  always  be  viewed  as  a  variation  of  parameters  as  in  the  above  example.  In  fact,  there  are 
examples  where  the  above  model  does  not  apply. 

Consider  the  pendulum  once  again.  Suppose  we  oscillate  the  pin  o'  the  pendulum  horizontally  instead  of 
vertically,  producing  motions  that  are  described  by 


xi  =  x2, 

x2  =  Csin(x1)  -  Bi2  +  aw2Dcos(xi)sin(wt). 


1 


(1.7) 

(1.8) 


Instead  of  the  moment  arm  sin(xi)  we  now  have  a  moment  arm  cos(xi)  and  the  naturally  occurring  feedback  is 
xv2Dcos(xx).  Linearization  of  this  system  of  equations  yields 


*1  ' 

'  0 

i 

Xi 

0 

.  *2  . 

c 

-B 

.  . 

+ 

au>2Dsm(wt) 

which  cannot  be  written  in  the  form  of  (1.6).  Consequently,  we  cannot  view  the  above  case  as  a  variation  of 
parameters. 

The  above  example  demonstrates  that  vibrating  a  system  component  does  not  always  produce  ” variation 
of  parameters’  as  in  the  vertically  vibrated  pendulum.  Consequently,  we  adopt  a  more  general  approach  that 
permits  the  analysis  of  problems  where  a  vibrated  system  component  may  result  in  nonlinear  functions  in  the 
governing  equations.  Consider  a  nonlinear  system 


i  =  /(*),  (1.10) 

with  an  equilibrium  point  at  the  origin  (  i.e.,  /(0)  =  0  ).  Vibrational  control  is  applied  by  oscillating  a  system 
component  or  process  at  high  frequency  and  low  amplitude.  For  instance,  in  the  case  of  a  jet  engine,  the  air- 

throttle  or  amount  of  fuel  injected  might  be  vibrated.  Let  h(wt)  =  sm(vjt)  denote  the  applied  high  frequency 

vibration.  It  is  assumed  that  the  vibration  affects  the  system  f(x)  through  some  naturally  occurring  feedback 
function  g(xy  w,  a),  which  depends  on  the  vibrated  component.  The  vibrationally  controlled  system  is  described 
by 

*  =  f(x)  +  h  {w^gix.w.a).  (1.11) 

For  convenience,  the  amplitude  of  h(xvt)  is  taken  to  equal  unity  and  the  amplitude  of  the  applied  vibration  is 
accounted  for  by  g(x,  w ,  a).  In  the  case  of  the  pendulum, 

/(x)  =  [x2,Csin(x!)  -  Bx2)t  (1.12) 

and  g(x,w,a)  =  [0,  aw2  D  sin(x  X)]T  for  the  vertically  vibrated  pin,  or  g{xyw,a)  =  [0,  a^2j5cos(xi)]r  for  the 
horizontally  vibrated  pin.  We  emphasize  once  again  that  g(x,w,a)  occurs  naturally,  and  is  not  measured  or 
computed  but  is  a  result  of  the  interaction  between  the  system  and  vibrated  component.  Obviously,  an  oscillating 
fuel  injection  rate  is  not  going  to  affect  the  jet  engine  in  the  same  fashion  as  an  oscillating  throttle.  Consequently, 
each  actuation  will  be  described  by  a  different  function  g(x,  w,a).  Since  g(x,w,  a)  depends  on  properties  of  the 
system  (which  are  fixed)  and  the  vibrated  component,  we  can  only  control  the  choice  of  the  component  to  oscillate 
and  the  frequency  and  amplitude  of  the  vibration.  This  choice  determines  the  form  of  <?(x,u>,a),  and  since  in 
certain  cases  there  exist  no  g(x,w,a)  that  will  allow  vibrational  control,  such  control  is  not  always  feasible. 

We  now  turn  to  the  question  of  stability.  Suppose  the  equilibrium  point  x  =  0  of  (1.10)  is  unstable,  and  that 
there  exist  one  or  more  accessible  system  components  or  processes  that  can  be  vibrated,  each  associated  with  a 
function  g(x,  w,  a)  that  is  known.  The  objective  of  the  theory  presented  in  this  paper  is  to  determine  a  stability 
criterion  for  (1.11).  Consequently,  if  a  certain  g(xywy  a)  satisfies  the  derived  stability  criterion,  then  oscillation 
of  the  corresponding  system  component,  with  specific  frequency  w  and  amplitude  a,  will  alter  the  stability  of  the 
system  and  result  in  vibrational  control.  Therefore,  the  developed  criterion  will  determine  if  vibrational  control 
is  feasible  for  various  accessible  system  components  or  processes  in  a  given  system. 

Vibrational  control  has  found  various  applications,  including  lasers  [2]  and  particle  beams  [3].  Initial  work  on 
developing  a  general  theory  of  vibrational  control  was  carried  out  by  Meerkov  [4].  He  discussed  the  effect  of  vibra¬ 
tional  control  upon  stability,  transient  motion  and  response  of  the  controlled  system.  In  subsequent  publications, 
several  specific  nonlinear  problems  were  discussed  [5],  but  no  general  vibrational  control  was  proposed.  Such  a 
theory  was  outlined  by  Bellman,  Bentsman  and  Meerkov  [1],  who  presented  criteria  for  the  control  of  nonlinear 
systems  by  linear  vibrational  control.  Further  nonlinear  results  are  discussed  in  [6],  including  conditions  for  and 
choice  of  stabilizing  vibrations. 

To  discuss  the  results  derived  in  [1],  consider  (1.11)  and  assume  that  the  Jacobian  matrix  df( 0)/dx  =  m  of 
/(x)  in  (1.11)  has  a  positive  trace.  A  classic  theorem  in  linear  algebra  states  that  the  trace  of  a  matrix  equals  the 
sum  of  its  eigenvalues  (see  for  example  [7,  p.251]  ).  Consequently,  if  the  trace  is  positive,  then  at  least  one  of  the 
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eigenvalue8  mus1  haw  a  positive  real  part  and  the  equilibrium  point  is  unstable.  This  does  not  imply,  however 

Edition  forTteabnUygaUVe  eqUiHbriUm  ^  iS  *  necessary  but  not  a 

.  ,fnellm^’  ®ent5man  and  Meerkov  [1]  only  considered  linear  vibrational  control,  which  limited  the  analysis 
to  hn  ar  functions  p(x,o  w)  =  Afx  in  (1.11).  They  proved  that  if  the  Jacobian  /'( 0)  has  a  positive  trace 
and  p(*,a  «;)  is  linear,  then  vibrational  control  is  not  feasible,  indicating  that  no  matrix  M  can  stabilize  the 
system  (1.11).  In  this  paper  we  consider  a  more  general  case  of  vibrational  control  via  a  nonlinear,  slowly 
iia  }/^8iiSlT’  *  J"  °!her  W°rdS’  WC  consider  Unctions  whose  rate  of  change  with  respect  to  x  is  bounded  (i.e. 
I.5’/.  x'j  —  w.  v  ®  s  ow  *n  tf1^s  case>  vibrational  control  may  be  possible  even  if  the  trace  of  the  Jacobian 
Matrix  is  positive.  Specifically,  it  will  be  shown  that  there  exist  nonlinear  functions  g(x,a,w)  that  stabilize  the 
system  (1.11)  even  if  its  Jacobian  /'( 0)  has  a  positive  trace. 

ctoivrf  fTn?-  thlS  ^aper  is  that  nonlinearities  in  9(x,a,w)  may  not  be  negligible,  and  can  affect  the 
stability  of  (1.11).  This  result  is  of  practical  importance  for  the  following  reason.  In  engineering,  it  is  common 

practice  to  linearize  a  system  before  analyzing  its  stability.  However,  if  a  Unear  system  is  considered,  then 
bellman,  Bentsman  and  Meerkov’s  result  indicates  that  vibrational  control  is  not  feasible  when  the  Jacobian  has 
a  positive  trace  (note  that  positive  traces  occur  in  a  wide  variety  of  engineering  systems  e.g.,  Uquid  rockets  |8l) 
Most  engineering  systems  are,  however,  nonlinear  and  it  is  possible  that  nonlinearities  in  g(x,  a,w)  may  stabilize 
the  system  even  if  its  Jacobian  trace  is  positive.  This  implies  that  one  should  not  discount  vibrational  control  for 
systems  that  exhibit  a  positive  trace.  Instead,  one  should  investigate  the  nonlinear  functions  g(x,  a,  w)  associated 
with  vibrational  open  loop  control  to  determine  if  they  satisfy  the  stability  criteria  derived  in  this  paper.  We  also 
note  that  the  theory  presented  in  this  paper  agrees  almost  exactly  with  numerical  solutions  (see  section  3  1) 


2  General  Derivation 

Consider  once  again  the  nonlinear  system 

i  =  /(*)  +  h(wt)g(x,w,a),  (2.1) 

where  h(wt)  =  sin(tirt),  x  €  IRn  is  the  state  space  vector  and  x  =  0  is  an  equilibrium  point  of  (1.10),  which  is  not 
necessarily  an  equilibrium  point  of  the  forced  system  (2.1).  It  is  assumed  that  /(x)  is  three  times  continuously 
differentiable  and  g{x)  is  four  times  continuously  differentiable. 

We  will  show  that  the  nonautonomous  system  (2.1)  can  be  approximated  by  an  autonomous  system 

y  =  F(y)-  (2.2) 

This  approximation  means  that  there  exists  a  function  u(t,  y),  which  is  small  for  all  time,  such  that  x(t)  = 
y(f)+u(f,y).  Consequently,  if  Y(t)  is  asolution  of  (2.2)  and  A '(f)  is  a  solution  of  (2.1),  then  A'«)-y(f)  =  u(t  Y(t)) 
is  small  for  all  time  t.  Approximately,  Y(t)  corresponds  to  the  time  average  of  X(t),  and  it  describes  the  slow 

reSP^K°f.  he  SySnem’  WhUe  corresponds  to  the  small  amplitude  high  frequency  system  oscillations 

excited  by  the  small  amplitude,  high  frequency  control  input.  In  essence,  there  exist  two  time  scales.  A  fast 

10  thC  ^gH  freqUency  contro1  inPut  and  the  suiting  high  frequency  system  response 
u(‘.  x  (t))i  the  slow  time  scale  describing  the  time  averaged  system  response  Y(t).  Since  Y(t)  is  a  slow  or 
averaged  response  it  is  described  by  a  time  averaged  equation.  In  the  case  of  vibrational  contro],  the  control 
mput  coupletl  with  the  system  response  u(t,Y(t))  yields  a  non-zero  average  that  can  stabilize  the  system. 

We  will  use  the  following  notation.  Since  w  and  o  are  constant,  we  will  express  g(x,w,a)  as  g(x).  Also,  we 
define  the  Jacobian  matrix  J  =  /'(0)  =  df{0)/dx  and  let  su 

p(x)  =  /(x)-Jx,  (2,3) 

<t>(vjt)  =  —e2Jg(0) sin(urf)  -  eg(0) cos(urt),  (2,4) 

where  e  =  \/w,  and  p(x)  is  the  sum  of  all  terms  of  second  order  and  higher  in  the  Taylor  expansion  of  /(x) 
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around  i  =  0.  Furthermore,  we  introduce  the  constant  vector  b 


b  = 


f  Jp{4>{rvt))dt  - 


£>9'(0)Jg(0) 


where  T  =  2n/w  and  y'(0)  is  the  Jacobian  matrix  of  y(x),  and  the  constant  matrix  A , 

c2  d(g'(y)Jg(y))i 


A  = 


dy 


:(°) 


(2.5) 


(2.6) 


where  g'(y)  is  the  derivative  of  g(y)  evaluated  at  y,  and  d{g'(y)Jg(y)](0)/dy  denotes  the  derivative  of  9'(y)Jg{y) 
evaluated  at  zero.  Finally,  we  let 


C  —  ^  d"  £  Mi  +  6$b\/E  +  68^/e  -f-  6q  +  6q6i  88q  65q6\  +  £^6q 

and  denote  a  ball  of  radius  6  centered  at  a  as  B(a,6). 


(2.7) 


Theorem  2.1  Consider  the  nonlinear  system  (2.1)  and  suppose  that  /( 0)  =  0,  ||y(0)||  <  w6q,  and  ||y'(£)||  <  w6\ 
for  all  £  £  £>(0,6).  Then  for  sufficiently  small  6,  6q  and  6\,  and  sufficiently  large  w,  there  exists  a  function 
u(t,y)  that  satisfies  the  following  properties:  ||u(f,y)||  <  2(6o  +  Mi)  for  all  t  and  for  all  y  £  B(0,6),  it  is  2i r/w 
periodic  in  t,  and  for  any  y  has  zero  mean  value.  Furthermore,  for  x(t)  governed  by  the  differential  equation 
(2.1),  y(t)  =  x(t)  -  u(t,y)  is  governed  by 

y  —  Ay  +  fc  +  O(C)  (2.8) 

for  all  y  £  B(0,6)  and  b ,  A  and  (  defined  in  (2.5),  (2.6)  and  (2.7),  respectively . 


While  a  detailed  proof  of  Theorem  (2.1)  is  given  in  the  Appendix,  an  outline  of  the  proof  is  provided  below.  A 
transformation  u(s,y)  is  constructed  that  satisfies  the  properties  of  the  theorem.  We  then  substitute  the  equation 
y(t)  =  x(t)  —  u(f,y)  into  (2.1)  and  bound  various  terms  so  that  we  can  rewrite  (2.1)  as  the  approximate  system 
y  =  F{t,y).  Next,  we  apply  the  method  of  averaging  to  derive  the  averaged  equation  y  =  Fov(y).  Linearization 
of  y  =  FQV{y)  at  the  origin  yields  the  result  of  the  theorem. 

The  analysis  in  this  paper  includes  Taylor  terms  up  to  second  order  in  60  and  6\.  Consequently,  the  resulting 
error  C  is  of  third  order.  If  higher  accuracy  is  desired,  then  more  Taylor  terms  can  be  included  although  more 
stringent  smoothness  constraints  will  be  imposed  because  we  will  have  to  ensure  that  higher  order  derivatives  exist 
for  the  functions  f(x)  and  g(x).  We  note  that  for  the  examples  considered,  a  second  order  analysis  is  sufficient 
and  is  in  excellent  agreement  with  numerical  integration  results  (see  example  3.1). 


2.1  Example:  The  Inverted  Pendulum 

Consider  the  vertically  vibrated  pendulum  described  by  (1.1)  and  (1.2).  These  equations  are  of  the  form  of  (2.1). 
Since  g(0)  =  0,  (2.4)  implies  that  <j>{wt)  =  0  and  (2.3)  show's  that  p( 0)  =  /( 0)  -0  =  0.  Consequently,  the  vector 
b  defined  in  (2.5)  equals  zero.  The  matrix  A  is  defined  in  (2.6)  and  can  be  expressed  in  the  following  form, 


g2  d(g'(y)Jgfa))„J  _  ,  g2  d  ( 

2  &y  WJ  2dy{ 


0 

aw2  D  cos(yi) 


'  0  1 

C  -B 

aw2  D  sin(yi) 


]} 


’  0  1 

d 

0 

0  1 

C  -  B 

&y 

a2w2D2  cos(yi)sin(yi)/2  _ 

C  -  (awD)2/ 2  -B 

Consequently,  Theorem  (2.1)  indicates  that  the  averaged  behaviour  of  the  system  is  governed  by 


in 

0  1 

Vi 

.  V2  . 

C-{awD)2/ 2  -B 

.  V2  . 

(2.10) 
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which  is  in  agreement  with  the  result  of  [4].  Note  that  the  term  C  —  (awD)7 / 2  is  negative  for  sufficiently  large  a 
or  u\  indicating  that  the  equilibrium  point  is  asymptotically  stable.  We  also  note  that  even  though  the  method 
in  this  paper  is  restricted  to  slowly  varying  g(x),  (i.e.,  ||f?'(x)||  <  w6\  <  w)  the  above  result  is  also  valid  for 
llp'OOII  it  w'  impose  the  slowly  varying  restriction  to  permit  inverting  the  matrbe  \I  +  uv]  in  (5.21).  In  the 
case  of  the  pendulum,  we  can  show  that  the  matrix  [/  +  uv]  has  an  inverse  even  if  ||(?'(x)||  it  w,  which  eliminates 
the  slowly  varying  restriction. 


3  Discussion  of  Results 

Theorem  (2.1)  implies  that  vibrational  control  can  result  in  an  equilibrium  shift.  For  such  a  shift  to  occur,  the 
vector  fc  defined  in  (2.5)  has  to  be  nonzero.  Equations  (2. 3), (2.4)  and  (2.5)  imply  that  such  an  equilibrium  shift 
can  occur  only  if  y(0)  is  nonzero.  In  this  case  there  are  two  possibilities.  The  first  possibility  is  that  the  average  of 
p(6(u>f))  is  nonzero.  Since  p(x)  is  defined  in  (2.3)  as  the  nonlinear  terms  of  /(x),  this  implies  that  nonlinearities 
in  f(x)  can  cause  an  equilibrium  shift.  Such  an  equilibrium  shift  would  be  of  order  O(||0||2)  =  0(<5g).  The  second 
possibility  is  that  the  term  p'(0)  Jg(0)  is  nonzero,  indicating  that  the  naturally  occurring  feedback  function  g(x)  can 
also  cause  an  equilibrium  shift.  In  this  case,  the  equilibrium  shift  would  be  of  order  0(e2  ||p'(0)  Jp(0)||)  =  O(Mi)- 
In  either  case,  if  the  equilibrium  shift  is  larger  than  6  our  analysis  fails  because  we  are  forced  outside  the  ball 

B(tU). 

Theorem  (2.1)  also  yields  a  useful  linear  result.  Consider  a  linear  system  of  the  form 

x  =  [J  +  sin(urt)£]x,  (3.1) 

where  \\B\\  <  w6\.  In  this  case,  g(x)  =  Bx  and  g'(x)  =  B.  Therefore  g(0)  =  0  and  we  can  set  So  ~  0  with  no  loss 
of  generality.  Application  of  Theorem  (2.1)  yields  the  averaged  equation 

V  =  [J7  e2BJB/2)y  +  O  ( 6(6  +  ^  +  S4/e)j  (3.2) 

However,  the  most  interesting  implication  of  Theorem  (2.1)  is  the  following:  the  operator  g’{y)Jg{y )  in  (2.6) 
is  a  nonlinear  operator  in  g(y).  Consequently,  nonlinearities  in  g(y)  may  result  in  linear  terms  in  the  averaged 
equation  (2.8),  and  can  influence  local  stability.  This  indicates  that  the  local  stability  of  the  nonlinear  system 
(1.11)  is  not  the  same  as  the  stability  of  a  corresponding  linearized  system.  It  is  possible  to  show  that  the 
nonlinearities  in  g(y)  can  alter  the  stability  of  a  system  with  a  positive  Jacobian  trace.  Stabilization  of  a  system 
with  a  positive  trace  is  illustrated  in  the  next  example. 


3.1  Example:  A  System  with  a  Positive  Jacobian  Trace 

In  this  example  we  consider  a  second-order  system  with  a  positive  trace.  Specifically,  we  consider  the  second-order 
system  derived  in  [8]  for  the  flow  potential  of  a  liquid  rocket  combustor, 


x  +  Ai±  +  A0x  =  0,  (3.3) 

where  x  is  a  nondimensional  flow  potential  perturbation  and  t  is  a  normalized  time.  In  an  unstable  liquid  rocket, 
unsteady  combustion  provides  negative  damping  that  drives  the  instability.  Since  the  damping  is  determined  by 
Ai,  negative  damping  corresponds  to  a  negative  coefficient  A\.  To  illustrate  the  point,  we  let  A\  -  -0.2  and 
Aq  =  1,  and  rewrite  (3.3)  as  the  following  second  order  system 


*1 

0 

1 

Xi  " 

0 

l 

Xi  ' 

.  *2  . 

i 

O 

-Ai  _ 

.  *2  . 

-l 

0.2 

.  *2  _ 

(3.4) 


The  Jacobian  matrix  of  (3.4)  has  a  positive  trace,  indicating  that  the  equilibrium  point  x  =  0  is  unstable. 

Bellman,  Bentsman  and  Meerkov  [1]  prove  that  it  is  not  possible  to  vibrationally  control  a  system  with  a 
positive  trace  if  the  function  g(x)  is  linear.  Consequently,  postulate  the  existence  of  a  nonlinear  function  g(x) 


9{x)  = 


0 

Q  +  0X 1X2 


(3.5) 


5 


that  describes  the  effect  produced  by  forcing  a  system  component.  We  stress  once  again  that  such  a  g{x)  would 
have  to  occur  naturally.  We  will  now  show  that  if  such  a  nonlinear  g{x)  exists,  it  will  stabilize  the  system  (we  do 
not  claim  that  such  a  g{x)  is  possible  in  rocket  motors).  A  discussion  of  the  reasoning  for  choosing  the  specific 
nonlinear  g{x)  given  in  (3.5)  is  provided  in  Section  (5.3)  in  the  Appendix. 

Given  the  above  choice  of  p(x),  we  write  the  forced  equation  as 


X\ 

0 

1 

X\ 

0 

.  . 

-Aq 

-A,  _ 

+ 

_  a  +  0xxx2 

sin(u;/). 


Let  w  —  70,  o  —  15  and  0  =  200.  For  these  values,  (3.6)  becomes 


X\ 

0  1 

Xj 

0 

.  . 

-1  0.2  _ 

.  *2  . 

+ 

15  +  200xjX2 

sin(urt). 


(3.6) 


(3.7) 


It  follows  from  (3.5)  that  ||p(0)||  <  a  <  w60  for  60  *  0.22.  Similarly  ||p'(x)||  <  0(xl  +  x2)  <  2 06  <  w6x  for 
6i  (5.72)6.  Consequently,  both  6q  and  6X  are  sufficiently  small  and  we  can  apply  Theorem  (2.1). 

We  need  to  calculate  the  vector  b  and  the  matrix  A  defined  in  (2.5)  and  (2.6),  respectively.  Notice  that  p(x), 
defined  in  (2.3),  is  zero  because  the  system  (3.4)  is  linear.  Consequently,  (2.5)  yields 


6=  - 


eV(0)Jp(0) 


However, 


<?'(*)  = 


0  0 
0x2  0x\ 


indicating  that  p'(0)  =  0,  which  implies  b  =  0  and  that  there  is  no  equilibrium  shift. 
Equation  (2.6)  yields  the  matrix  A , 


J- 


f2  d{g'(y)Jg{y)) 
2  dy 


=  J-H 


(0) 

0 


0 

0 

r  o  i  1 

0 

\ 

J  2  dy{ 

.  Pyi 

Pyx  . 

— i 

1 

o 

1 

- 1 

q  +  0yxy2  _ 

J 

2  dy  [  aPyi  -  Axa0yi  +  02y\y2  -  Ax02y\y2  J  [  -Ao  +  e2AXQ0/2  -A\  -  e2a0/2  J  ‘ 
Consequently,  Theorem  (2.1)  implies  that  the  averaged  motion  of  the  system  is  governed  by 


Vi  * 

0  1 

'  2/1  * 

2/2 

-Aq  +  e2Ax  q/3/2  -Ax-e2a0/2 

-  2/2 

Substituting  the  numerical  values  for  Aq,  Ax,  a,  0  and  e  -  l/w  yields 


(3.8) 

(3.9) 


(3.10) 


(3.11) 


yi 
.  V2 

which  is  asymptotically  stable. 

Since  the  solution  X(t)  of  (3.6)  is  given  by  X(t)  =  Y(t)  +  u(t,Y{t )),  where  Y{t)  is  a  solution  of  (3.12)  and 
tends  towards  the  origin  as  time  tends  to  infinity,  X{t)  must  remain  close  to  the  origin  for  all  time  because 
u(t,Y(t))  is  small  for  all  time.  The  construction  of  u(t,y),  as  defined  in  the  Appendix  (see  (5.4), (5.8)  and  (5.9)), 
implies  that  if  $(0)  ^  0  then  u(t,y)  -/+  0  as  y  — ►  0.  In  this  case  p(0)  ±  0,  indicating  that  u{t,Y(t))  does  not 
converge  to  zero  as  y(t)  tends  to  zero.  Consequently,  X(t)  remains  close  to  zero  for  all  time  but  does  not  tend 
to  zero  as  time  goes  to  infinity.  Strictly  speaking,  the  equilibrium  point  x  =  0  of  (2.1)  is  not  asymptotically 
stable;  indeed  x  =  0  is  not  an  equilibrium  point  but  is  the  center  of  a  small  asymptotically  stable  limit  cycle. 
This  limit  cycle  is  the  asymptotically  stable  orbit  X(t)  =  u(t,0)  ^  0.  We  refer  to  x  =  0  as  a  slow  equilibrium 
point  because  y  =  0  is  an  equilibrium  point  of  the  slow  or  time  averaged  system  (2.8),  and  we  say  that  x  =  0  is 


0  1 
-1.061  -0.106 


Vi 

.  V2 

(3.12) 


6 


slowly  asymptotically  stable  because  the  equilibrium  point  y  =  0  of  the  slow  system  (2.8)  is  asymptotically  stable. 
When  we  refer  to  slow  equilibrium  points  or  slow  stability,  we  refer  to  the  properties  of  the  time  averaged  system 
(2.8).  The  true  dynamics  are  small  oscillations  about  the  slow  or  averaged  dynamics  and  hence  display  the  same 
qualitative  behaviour.  Prom  a  practical  point  of  view  we  have  achieved  our  control  objective  to  keep  the  system 
(2.1)  in  a  small  neighbourhood  of  the  origin.  Therefore,  if  there  exists  an  accessible  component  in  a  liquid  rocket 
motor  that  can  produce  a  naturally  occurring  feedback  function  g(x)  =  (07a  +  /Jx^]7,  then  we  can  achieve 
vibrational  control  by  vibrating  this  component. 

It  is  interesting  and  instructive  to  compare  results  obtained  by  this  analysis  with  a  numerical  simulation.  We 
can  analytically  solve  the  time  averaged  equation  (3.12)  to  derive  the  following  analytic  expression  for  Vi(<) 

yi(0  =  c-ao5^  (3.13) 

where  Y\(0)  is  the  initial  displacement  and  V^O)  is  the  initial  velocity.  Figure  (1)  compares  Ki(t)  of  (3.13)  with 
a  X\(t)  calculated  by  numerically  solving  (3.7).  Since  the  initial  conditions  for  the  slow  solution  Y(t)  are  not 
known,  they  are  matched  to  the  initial  conditions  shown  by  the  numerical  simulation.  Figure  (1)  shows  that 
the  slow  equilibrium  point  x  =  0  of  the  forced  system  (3.6)  is  indeed  slowly  asymptotically  stable  (i.e.,  X\(t) 
approaches  a  small  asymptotically  stable  limit  cycle)  but  is  not  asymptotically  stable  (A'i {t)  y*  0).  Furthermore, 
Fig.  (1)  shows  excellent  agreement  between  the  behaviour  predicted  by  the  developed  theory  and  the  numerical 
simulation. 

4  Conclusion 

In  this  paper  wre  present  a  criterion  for  nonlinear  vibrational  open  loop  control.  Previous  work  that  was  restricted 
to  linear  control  is  extended  to  include  analysis  of  nonlinear,  vibrational  control.  It  has  been  previously  shown 
that  linear  vibrational  control  is  not  feasible  if  the  Jacobian  matrix  has  a  positive  trace.  This  paper  demonstrates 
that  nonlinear  vibrational  control  is  possible  even  if  the  trace  of  the  Jacobian  is  positive.  This  result  is  significant 
because  a  large  number  of  nonlinear  engineering  systems  exhibit  a  positive  Jacobian  trace  and  yet  may  be  stabilized 
by  nonlinear,  open  loop,  vibrational  control.  Finally,  it  is  shown  that  the  theory  developed  in  this  paper  is  in 
excellent  agreement  with  numerical  results. 


5  Appendix 

In  this  section  we  prove  Theorem  (2.1)  and  discuss  the  corresponding  change  of  variables  x(t)  =  y(t)  +  u(t,y). 
We  begin  by  assuming  that  the  investigated  system  is  described  by 


x  =  f(x)  + h(wt)g(x,w,a),  (5.1) 

where  x(t)  €  IRn,  /  €  C3(Rn,Rn),  /( 0)  =  0,  h(wt)  =  sin (wt),  w  »  1,  and  g  €  C4(Rn  x  R  x  R,R").  We  perform 
a  local  analysis  that  will  be  restricted  to  a  ball  of  radius  6  centered  at  the  origin.  In  addition,  since  xv  and  a  are 
constant,  we  write  g(x,  w,  a)  simply  as  g(x)  and  impose  the  following  smoothness  constraints 


ll/'(0)|| 

< 

p, 

0 

< 

O, 

ns(o)ii 

< 

w60, 

0 

< 

<$0, 

wm 

< 

w6u 

0 

< 

where  /'(x)  denotes  the  derivative  of  /  evaluated  at  x  and  f  €  B(0,  <5).  To  simplify  the  algebra,  we  introduce  a 
fast  time  variable  a  =  wt,  define  e  =  \/w,  denote  dx/ds  as  x  and  rewrite  (5.1)  in  the  fast  time  scale 


x  =  ef(x)  +  eh(s)g(x). 


(5.3) 
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5.1  The  Transformation 


To  prove  Theorem  (2.1)  we  introduce  the  change  of  variables  i(s)  =  y(s)  +  u(s,y).  Next,  we  define  the  function 
u(s,y)  and  determine  some  of  its  properties;  that  is, 


u(s,  y)  =  a{y)  sin(s)  +  (3{y)  cos(s), 


(5.4) 


where  a(y),/3(y)  €  lRn.  The  functions  a(y)  and  0{y)  are  chosen  so  that  u(s,y)  satisfies  the  partial  differential 
equation 

«a(s,y)  =  eJu(s,y)  +  eh(s)g{y),  (5.5) 

where  J  -  /'( 0)  is  the  Jacobian  matrix  and  the  subscript  s  denotes  a  partial  derivative  with  respect  to  s.  Note 
that  for  any  fixed  y  the  above  equation  is  an  ordinary  differential  equation  in  u(-,y).  Substituting  (5.4)  into  (5.5) 
and  equating  the  coefficients  of  the  sines  and  cosines  yields 


-  P{y)  -  eJ<*{y)  =  eg(y), 

(5.6) 

a(y)  -  eJ0(y)  =  0. 

(5.7) 

0(y)  yields 

Q(y)  =  -c2[/  +  c2J2]-1Jy(y), 

(5.8) 

0(y)  =  —eJa(y)  -  eg{y), 

(5.9) 

where  the  inverse  matrix  [7 -he2  J2]  1  is  well  defined  provided  e  is  small  enough  to  satisfy  the  inequality  ||e2  J2||  <  1. 
To  derive  approximate  equations  for  a(y)  and  0(y)  we  need  the  following  bound  on  g{y), 


llff(y)ll  <  ||ff(0)  +  g{y)  -  s(0)||  <  ||y(0)||  +  ^||y||  < 


(5.10) 


which  holds  for  all  y  €  B{0,6).  Next,  we  represent  the  inverse  matrix  [I  +  e2J2]  1  as  the  geometric  series 

[/  +  r2  J2]"1  =1-  e2J2  +  e4J4  -...  =  /  +  O(eV).  (5.11) 

Using  (5.8)  through  (5.11)  yields  the  following  approximate  expressions 


a(y)  =  -e2Jg{y)  +  O(e360  +  e3661), 

0(y)  =  -c9(y)  +  c3J29(y)  +  0(c4«o  +  f4Wi). 


(5.12) 

(5.13) 


To  complete  the  discussion  of  the  properties  of  u(s,y),  we  need  bounds  on  u(s,y)  and  the  partial  derivative 


uv(s,y).  We  begin  by  bounding  the  inverse  matrix  [7  +  e2J2]  l.  Equation  (5.11)  implies 


II  [I  +  c*J 


2  t2  1  —  1 1 


+  +  . . .  <  1  +  eV  +  . . .  <  1/(1  -  eV4). 


^2  _2\ 


It  follows  from  (5.8),  (5.9),  (5.10)  and  (5.14)  that 

My)ll  < 

mv)\\  < 


ea(6o  +  (56i) 

1  —  e2c72  ’ 

6q  +  66, 

1  -  e2a2 


(5.14) 

(5.15) 

(5.16) 


for  all  y  €  73(0, 6).  To  derive  the  desired  bound  on  u(s,  y)  we  need  only  to  note  that  (5.4)  implies  ||u||  <  ||q||  + 1|/?||, 
indicating  that 

||«(«,y)||  <  +  <  2 (So  +  66,)  =  O(60  +  66,), 

I  — * 


(5.17) 
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which  holds  for  all  a  and  sufficiently  small  e.  The  bound  on  uv(s,y)  =  o'(y)sin(s)  +  0'{y)cos(s)  is  also  straight¬ 
forward.  Since  (5.8)  and  (5.9)  imply 

e'(y)  =  -e2{I  +  e2J2}-1Jg'(y),  (5.18) 

0\y)  =  - eJa'(y)~  eg'(y ),  .  (5.19) 

then  using,  (5.10),  (5.14),  (5.18)  and  (5.19)  one  obtains 

K(«,»)ll<  =  Q(^i )»  (5.20) 

which  holds  for  all  s. 

5.2  Proof  of  Theorem  2.1. 

We  begin  by  noting  that  the  transformation  u(s,y)  constructed  in  the  previous  section  satisfies  the  constraints 
outlined  in  the  theorem.  The  transformation  x(s)  =  y(s)  +  u(s,  y)  implies  dx/ds  =  x  =  y  +  u,  +uvy.  Substituting 


this  relationship  into  (5.3)  yields 

[7  +  ttv(s,y)]y  +  tt,(s,y)  =  ef{y  +  tt)  +  sk{s)g{y  +  u).  (5.21) 

Equation  (5.20)  implies  ||uv(s,y)||  <  1  for  sufficiently  small  Su  for  all  y  e  B(0,6)  and  for  all  s.  Consequently,  the 
inverse  matrix  [I  +  uv(s,y)]-1  is  well  defined  and  we  can  rewrite  (5.21)  as 

y=[I  +  uv(s,y)}-l[ef(y  +  u)  +  eh{s)g{y  +  u)-us(s,y)}.  (5.22) 

The  following  relationships  will  be  used  to  simplify  (5.22), 

p(x)  =  /(x)  -  Jx,  (5.23) 

9(^»u)  =  9(y  +  u)-g(y)-g'(y)u.  (5.24) 

where  p(x)  is  defined  as  before  and  q(y,  u)  represents  the  sum  of  all  terms  of  second  order  and  higher  in  the  Taylor 
expansion  of  g(y  +  ti)  around  u  =  0.  It  follows  that 

P(0)  =  0,  p'(0)  =  0,  p(x)  =  0(||x|i2),  (5.25) 

q(y,  0)  =  0,  gu(y,  0)  =  0,  q(y,u)  =  0(||u||2).  (5.26) 

Using  (5.23)  and  (5.24)  we  can  rewrite  (5.22)  as 


y  =  l1  +  uv(s-  2/)]  1  \£jy  +  y)  +  ep(y  +  it)  +  eh{s)g{y)  +  eh(s)g'{y)u{s,  y)  +  eh(s)q{y ,  it)  -  it,  (a,  y)].  (5.27) 

Substituting  (5.5)  into  (5.27)  yields 

y  -  [I  +  uv(s>  p)]_1  [£  Jy  +  ep(y  +  *0  +  ^{s)g'(y)u{s,  y)  +  eh{s)q{y,  it)].  (5.28) 

Approximating  the  inverse  matrix  [7  +  uv(a,y)]_1  as  a  two  term  series  with  a  second  order  error, 

I/  +  “v(s,J/)r1  =I-Uy(s, y)  +  0(||uj2)  =7-  uv(s,y)  +  0(6?),  (5.29) 

and  substituting  (5.29)  into  (5.28)  yields 

y  =  s[7  -  xiy(s,  y)  +  0(^)][Jy  +  h(s)g'(y)u(s,  y)  +  p(y  +  it(a,  y))  +  h(s)q{y, tt(s,  y))]  =  eF(s,  y).  (5.30) 

We  are  now  in  a  position  to  apply  the  method  of  averaging.  Since  Fva,y)  is  periodic  in  a  with  a  period  2tt, 
we  can  approximate  the  non-autonomous  system  y  =  <rF(s,y)  as  the  autonomous  averaged  system  y  =  sFot,(y) 
where 

1  f 2* 

Fav(y)  =  —  y0  F(T,  2/)*,  (5.31) 
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(see  [9,  p-4 1 2]  for  a  discussion  of  averaging).  Consequently,  the  averaged  equation  is  given  by 
e  f 2* 

y=2^J0  l1 -uv(T<y)  +  0(62l)}\Jy  +  h(T)9'(yMT,y)  +  p{y+u{T,y))  +  h{T)q{y,u(T,y))}dT.  (5.32) 

Expanding  (5.32)  yields 

y  =  yr  jfl  Jy  +  h(T)9'{yMT,y)  +  p{y  +  u{r,y))  +  h(T)q{y,u{T,y))  -  uv(T,y)Jy 

~  ^y(T,y)h{r)g'(y)u(T,y)  -  uv(r,y)p(y  +  u(r,y))  -  uv(r,  y)h(r)q(y,  u(r,y)) 

+  O(66f  +  6o6f/c  +  66i/e  +  6$6f)  ]dr.  (5.33) 


The  terms  Uy(r,y)Jy  and  ^(r,  y)h(T)y'(y)u(T,y)  consist  of  an  odd  number  of  sinusoidal  functions  and  thus 
average  to  zero.  The  term  Jy  is  constant  with  respect  to  r  and  can  be  taken  outside  the  integral.  Finally,  since 
h(s)  =  sin(s)  and  u(s,y)  =  a(y)sin(s)  +  /3(y)  C06(s),  then  averaging  the  term  h(T)g'(y)u(T,y)  yields 


Jo  h{T)9'(y)u(T,  y)dr  =  ^  J  sin2(r)y'(y)a(y)  +  sin(r)  cos(t) g'(y) (3 (y) dr  =  g’{y)a{y)/ 2.  (5.34) 


Using  the  approximate  expression  (5.12)  for  a(y)  in  (5.34)  lets  us  rewrite  (5.33)  as 

,  s  f 

y  =  £Jy-—g\y)Jg{y)+—  [p(y+u(T,y))+/i(r)9(y,u(T,y))-uv(T,y)p(y+u(T,j/))-uv(r1y)h(r)9(y,u(r,y))]<iT 

+sO{e2606l  +  66]  +  60S  3Je  +  SSj/e  +  (5.35) 

To  complete  the  proof  we  have  to  bound  the  integral  in  the  (5.35).  The  bounds  on  uv(t,  y)p(y  +  tt(r,y))  and 
uv(r,y)h{T)q(y,u{T,y))  follow  from  (5.17),  (5.20),  (5.25)  and  (5.26);  that  is, 

wv(T>y)p(l'  +  u(T-l/))  =  °(IKIIIIy  +  u||2)  =  0(626i  +  6606:  +  (5.36) 

uv{T,y)h(T)q(y,u{T,y))  =  0(||u„||||u||2)  =  0(<5q6i  +  +  <52^).  (5.37) 


In  order  to  get  bounds  on  the  remaining  terms,  p(y  4-  u(r,  y))  and  h(r)q(y,  u(t,  y)),  we  will  require  the  following 
notation.  Denote  the  second  order  Taylor  expansion  of  g(y  -f  it)  at  u  =  0  as 


g(y  +  u)  =  g(y)  +  g'(y)u  +  ly"(y)(u,  u)  +  0(||u||3),  (5.38) 

where  (u,ti)  denotes  a  tensor  and  gn{y)  is  the  corresponding  three  dimensional  array  of  coefficients  evaluated  at 
y .  It  follows  that  q{y,u)  =  g"{y){u,u)/ 2  +  0(||u||3).  Consequently,  the  average  of  h(r)q{y,u(T,y))  is  written  as 

2tt  Jo  h(TWy' T ’  y^dr  =  ^  JQ  h(T)9”(y) («(T>  V)i  u(t,  y))dr  +  0(||u||3).  (5.39) 

Since  each  term  of  ^(T)(u(T,t/),u(r,  y))  consists  of  an  odd  number  of  sinusoidal  functions,  the  resulting  average 
is  zero.  Hence,  (5.39)  is  reduced  to 

h(r)q{y ,  u(r,  y))dr  =  0(6^  +  +  63<5?).  (5.40) 

With  the  aid  of  the  bounds  (5.36),  (5.37)  and  (5.40)  we  can  rewrite  (5.35)  as 
y  =  eJy-  ^g'{y)Jg{y)  +  ^  ^  p(y  +  U(T> y^dr  +  e0^2^\  +  +  M?/e  +  M\fe  +  +  ^i)-  (5'41 ) 
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Equation  (5.41)  is  of  the  form  y  -  F(y)  +  eO(. . .)  where 


£*3  - 

F(y)  =  eJy  -  -^g'iyjJgiy)  +  —  j  p{y  +  u{T,y))dT. 

Since  we  are  concerned  with  local  behaviour  at  the  origin,  we  linearize  (5.41)  about  y  =  0,  to  get 


y  =  F( 0)  +  —(0)  y  +  eO(62  +  ...). 


Expanding  the  above  yields 


d(g'(y)J9(y)), 


(0)  y  + 


g(^/oP(»  +  u(T.!/))^) 

dy 


(0)  v 


+eO(62  +  e2606l  +  66\  +  60b\/e  +  66\/c  +  «§  +  66^  +  %6i).  (5.44) 

We  now  complete  the  proof  by  bounding  the  last  term  in  (5.44).  Since  the  derivative  of  p(x)  exists  and  is 
continuous  by  assumption,  we  can  move  the  partial  derivative  d/dy  inside  the  integral  to  get 


5[p(y  +  u(r,j/))] 


(0)  }  dr 


where 


d\p(y  +  u(t,j/))| 


(0)  =  p'(u(r,  0))  +  p'(u(r,  0  ))uv(r,  0). 


Since  p(x)  €  lRn,  then  p'(a)  £  Rnxn  is  a  matrix  valued  function.  Letting  denote  the  ij’th  element  of  the 
matrix  M  and  Aj(a)  =  [p'(a)]>j  €  R,  then  using  (5.46)  lets  us  write  the  if  th  term  of  (5.45)  as 


9  (sr  fo*p(y  +  u(r-  i/))*) 

dy 


■  e  [2* 

“2 n  J0  ^^(T’0))  +  (^k(u(r,0))[tiv(T,0)]*j)d7 


where  the  tensor  notation  ( )  implies  a  summation  over  the  index  k.  Expanding  A,  and  A*  as  first  and  zero  order 
Taylor  series  about  the  origin  yields 


, 9 Ionp(y  +  u(T’ y))dr ) J  c  [2r , 

dy  ^  =  2rr  J0  +  4j(°Mt>°)  +  (MO^u^O)]*.,)}  dr+ffO(||u||2+||u||||uv||). 

*■  -1  v 

Equation  (5.25)  implies  Aj(0)  =  0  and  the  averages  of  ^(0)u(r,  0)  and  A* (0)K(r,  ())]*,•  are  zero.  Consequently! 

[5(^/o2,p(!/  +  «(T.I/))dr]  1 

I  _ 1 _ /  /A\  _ It-  II  tl  II 2  .  II  tin  rtn  nv  -  _  _  «  « 


(0)  y  =  eO(||y||||u||2  +  ||y||||ti||||itv||)  =  0(66$  +  A506i  +  626j).  (5.49) 


The  bound  (5.49)  allows  us  to  rewrite  (5.44)  as 


0  ))dr  — 


+  eJy 


f3  \d(g,(y)Jg(y))i 


(0)1  y-h 


tO(62  +  £2606i  +  662  +  6062/e  +  66\/e  +  «5q  +  6^1  +  +  A50<5i). 


;?  a.  X„X3  /«-  J_  £*4  le.  A  x3  .  elc  ,  c c2 
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According  to  the  definition  (2.4) 
then  (5.12)  and  (5.13)  imply 


<£(t)  =  —e2Jg(0)  sin(r )  -  ey(0)  cos(t), 
<P(t)  =  u(r,0)  +  0(e%- +  e3^). 


Using  (5.52)  in  (5.50)  yields 


y  = 


r7* 

/  P(0(T ))dr  — 

Jo 


rV(0)Jp(0) 


e*3  r 


+  eJy-- 


d(g'(y)Jg(y)) 


dy 


(0) 


y+ 


eO{62  +  e2606:  +  662  +  <5063/<r  +  66\/e  +  +  %6X  +  Mg  +  5<506i  +  e36g). 

Rewriting  (5.53)  in  the  original  time  scale  t  yields 


y  =  Ay  +  b  +  O(C). 

where  A,  b  and  £  are  as  defined  in  the  theorem.  □ 


(5.51) 

(5.52) 


(5.53) 

(5.54) 


5.3  Choice  of  g(x)  in  Positive  Trace  Example 

In  equation  (3.5)  we  let  y(x)  =  [0,a  +  0x\x2 J7-.  This  hypothetical  choice  of  g(x)  is  not  arbitrary.  We  know  that 
the  sign  of  A\  creates  an  instability.  Consequently  we  wish  to  change  the  sign  of  this  coefficient  by  applying 
vibrational  control.  Consider  the  averaged  equation  (2.8).  If  we  denote  the  vector  g'(x)Jg{x)  as  [Gj(x),  G2(x)]T 
then  the  matrix  A  defined  in  (2.6)  can  be  written  as 


A-j-kVMg'Wm.j.k 


dG j(0)/3xi  dGx{ 0)/dx2  ‘ 
8G2(  0)/dxx  dG2{  0)/8x2 


(5.55) 


where  k  is  a  positive  constant.  For  A  to  have  a  negative  trace  either  dGi(0)/dxx  or  5G2( 0)/dx2  must  be  positive, 
or  both.  Consequently,  letting  dG2{0)/dx2  =  c  be  a  positive  quantity  implies  G2(x)  =  cx2.  It  follows  that 

G2(x)  =  ~~  +  0-2<?2(x)^~-  =  cx2.  (5.56) 

If  we  consider  the  first  term  only,  we  can  set 


f  ,dg2 

P2(I)_  =  CI2. 


(5.57) 


Equation  (5.57)  is  a  partial  differential  equation  in  g2(x),  which  can  be  solved  by  separation  of  variables.  Unfor¬ 
tunately,  the  solution  to  (5.57)  is  g2(x)  =  ici-^/jxixTI  which  is  singular  at  the  origin  and  violates  the  assumption 
that  g(x)  is  continuously  differentiable.  Consequently,  we  let  g2(x)  =  a  +  fix \x2,  approximating  the  square  root 
dependance  of  g2(x)  near  the  origin.  If  we  now  set  g^x)  =  0,  then  g(x)  =  [0,a  +  f3x1x2}T.  It  is  noteworthy 
that  the  last  term  in  (5.56)  suggests  that  g2(x)  =  Kx2  might  also  be  a  viable  feedback  function.  Such  a  choice 
requires,  however,  that  K  >  w,  which  violates  the  assumption  that  ||y'(x)||  <  wSi  <  w. 
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